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The  Davis-Stenger  methodology  was  adapted  to  the  problem  of 
vibration  control  of  flexible  structures.  The  spectral  factorization 
methodology  avoids  the  difficult  numerical  problems  associated  with 
the  solution  of  the  Riccati  partial  differential  equations  which  arise 
in  the  time  domain  approach  for  designing  stabilizing  controllers.  In 
this  way  distributed  phenomena,  like  travelling  waves,  which 
characterize  the  macroscopic  dynamics  of  flexible  structures  are 
retained  in  the  model,  and  their  interaction  with  the  control  system 
is  preserved  in  the  analytical  design  process.  lomputational 
algoritnms  were  developed  and  several  prctotj'pe  systems  were  treated 
including  the  Euler  Beam  and  a  simple  two  dimensional  system. 


Second  part  of  the  research  involved  t^  use  of  a  mathematical 
technique  for  asymptotic  analysis  called  '"^mogenizatioi^'^ originally 
developed  by  I.  Eabuska,  to  produce  simplified  models  Qfor  flexible 
structures  with  a  regular  (periodic)  infrastructure,  '^'omogenization 
cf  the  model  for  a  structure  with  a  regular  infrastructure  produces  a 
model  with  smoothly  varying  "Effective  parameters  for  mass  density, 
local  tendon,  and  damping  that  represents  a  flexible  structure  with  a  ^  y 
uniform  ""homogenised'‘*^nternal  structure. 

The  homogenization  technique  does  not  require  the  a  priori 
assumption  of  a  specific  continuum  structure  as  the  approximation  for 
a  given  lattice  structure.  Instead,  the  asymptotic  analysis  of  the 
original  structure  produces  the  distributed  continuum  approximation 
model  of  the  lattice  structure  in  the  limit  as  some  characteristic 
parameter  (e.g.,  inter-cell  dimension)  in  the  structural  model  gees  to 
zero.  Moreover,  the  natural  averaging  process  is  developed  in  the 
course  of  the  analysis.  It  is  easy  to  construct  examples  in  which  the 
procedure  of  averaging  parameters  over  a  c' eracteristic  volume  leads 
to  incorrect  approximations  for  the  system  dj.tcamics.  The 
homogenization  methods  used  in  this  re.  earoh  are  based  on  the 
assumption  of  a  periodic  infrastructure  in  the  original  model.  This 
is  net  necessary,  and  random  structures  can  also  be  treated,  if  the 
randomness  has  sufficient  ergodicity  properties  (in  the  spatial 
variables). 

Homogenization  and  asymptotic  analysis  can  also  be  carried  out  in 
the  context  of  control  and  state  estimation  problems  for  heterogeneous 
structures.  It  is  important  that  the  control  and  homogenization 
procedure  not  be  done  separately,  since  one  can  construct  examples  in 
which  control  designs  based  on  averaged  models  are  not  correct  as 
approximations  to  the  optimal  (e.g.,  regulator)  control  laws  for  the 
original  problem.  'rfhile  control  and  filtering  theory  with 
homogenization  is  not  very  advanced  at  this  stage,  prototype  problems 
can  be  analyzed  to  a  point  where  the  basic  features  of  design 
algorithms  are  clear. 
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Suomary  of  Phase  I  Research 

Interaction  of  the  control  system  with  the  structural  dynamics  of 
the  physical  system  is  one  of  the  fundamental  Issues  in  large  space 
structure  applications.  Our  work  is  intended  to  contribute  directly 
to  understanding  this  interaction  by  using  models  which  capture  the 
essential  distributed  character  of  the  system,  and  using  analytical 
techniques  which  preserve  the  character  of  the  physical  system  in  the 
model  simplification  process.  The  methods  we  have  used 
Yiener-Hopf/spectrml  factorization  methods  for  design  of  distributed 
control  systems  and  homogenization/asymptotic  analysis  for  model 
simplification  —  have  tremendous  potential  for  the  analytical 
treatment  of  complex  structural  control  problems,  including  the 
sysnthesls  of  computer-aided-design  methods  for  large  space 
structures.  In  Phase  I  of  this  project,  we  have  concentrated  on  the 
treatment  of  a  few  simple  prototype  systems.  Further  work  is  needed 
to  adapt  and  enhance  the  methods  to  treat  complex  structures.  The 
analytical  methods  themselves  do  not  require  substantial  extensions. 
Rather,  their  potential  for  the  treatment  of  complex  flexible 
structures  should  be  developed. 

The  main  emphasis  in  the  first  phase  of  this  work  has  been  the 
adaptation  and  enhancement  of  certain  Wiener-Hopf  methods  for  control 
system  design  used  by  J.  Davis  for  the  treatment  of  linear,  dynamic, 
distributed  parameter  models  of  flexible  structures.  Davis  developed 
a  frequency  domain  methodology  for  computing  optimal  (regulator) 
feedback  gains  for  linear  distributed  parameter  control  systems  by 
Viener-Hopf  spectral  factorization.  The  numerical  algorithms  for 
executing  the  spectral  factorization  were  based  on  some  earlier  work 
of  F.  Stenger.  Ve  have  adapted  the  Davis-Stenger  methodology  to  the 
problem  of  vibration  control  of  flexible  structures.  A  generic 
problem  of  this  type  is  the  figure  control  of  a  large  space  antenna. 
Ve  have  carried  out  the  analysis  and  computed  the  optimal  feedback 
regulator  control  laws  for  several  examples  including  the 
Euler-Bemoulli  beam  model  and  a  two-dimensional  prototype 
(experimental)  system  studied  by  J.  Lang  and  D.  Staelln. 

This  portion  of  the  research  has  demonstrated  the  effectiveness 
of  frequency  domain  —  spectral  factorization  methods  for  the  design 
of  control  and  state  estimation  algorithms  for  flexible  structures 
described  by  linear  distributed  parameter  models  (hyperbolic  partial 
differential  equations).  In  this  approach  it  is  not  necessary  to 

reduce  the  models  to  finite  dimensional  (lumped  parameter)  models  at 
the  outset  of  the  design  procedure.  The  infinite  dimensional 

character  of  the  system  is  preserved  throughout  the  design  process. 
The  spectral  factorization  methodology  avoids  the  difficult  numerical 
problems  associated  with  the  solution  of  the  Rlccati  partial 

differential  equations  which  arise  in  the  time  domain  approach  for 

designing  stabilizing  controllers.  In  this  way  distributed  phenomena, 
like  travelling  waves,  which  characterize  the  macroscopic  dynamics  of 
flexible  structures  are  retained  in  the  model,  and  their  Interaction 
with  the  control  system  is  preserved  in  the  analytical  design  process. 


.%  s 


In  the  second  part  of  the  research  ve  have  examined  the  use  of  a 
mathematical  technique  for  asymptotic  analysis  called 
"homogenization",  originally  developed  by  I.  Babuska,  to  produce 
simplified  models  for  flexible  structures  with  a  regular  (periodic) 
Infrastructure.  Homogenization  of  the  model  for  a  structure  with  a 
regular  Infrastructure  produces  a  model  with  smoothly  varying 
"effective"  parameters  for  mass  density,  local  tension,  and  damping 
that  represents  a  flexible  structure  with  a  uniform  "homogenized" 
internal  structure.  The  derivation  of  continuum  models  for  complex 
structures  with  a  regular  infrastiructure  has  been  studied  for  many 
years  in  applied  mechanics.  In  most  cases  the  continuum  models  are 
baaed  on  local  averages  of  the  physical  parameters  (e.g.,  mass 
density,  stress,  strain,  etc.)  over  some  characteristic  volume  of  the 
structure.  The  averaged  quantities  computed  in  this  way  are  related 
to  the  associated  quantities  in  a  postulated  continuum  structure.  For 
example,  the  mass  density  and  stress  tensors  in  a  long  truss  with  a 
regular  lattice  structure  have  been  related  to  the  distributed 
parameters  in  a  beam  (in  the  work  of  Noor,  Nayfeh,  and  Renton,  among 
others). 

Our  technique  does  not  require  the  a  priori  assumption  of  a 
specific  continuum  structure  as  the  approximation  for  a  given  lattice 
structure.  Instead,  the  asymptotic  analysis  of  the  original  structure 
produces  the  distributed  continuum  'approximation  model  of  the  lattice 
structure  in  the  limit  as  some  characteristic  parameter  (e.g., 
inter-cell  dimension)  in  the  structural  model  goes  to  zero.  Moreover, 
the  natural  averaging  process  is  developed  in  the  course  of  the 
analysis.  It  is  easy  to  construct  examples  in  which  the  usual 
procedxire  of  averaging  parameters  over  a  characteristic  volume  leads 
to  Incorrect  approximations  for  the  system  dynamics.  The 
homogenization  methods  used  in  our  research  are  based  on  the 
assumption  of  a  periodic  infrastructure  in  the  original  model.  This 
is  not  necessary,  and  random  structures  can  also  be  treated,  if  the 
randomness  has  sufficient  ergodlcity  properties  (in  the  spatial 
variables).  Numerical  evaluations  of  the  averaged  model  are  more 
difficult  in  this  case. 


Homogenization  and  asymptotic  analysis  can  also  be  carried  out  in 
the  context  of  control  and  state  estimation  problems  for  heterogeneous 
structures.  It  is  Important  that  the  control  and  homogenization 
procedure  not  be  done  separately,  since  one  can  construct  examples  in 
which  control  designs  based  on  averaged  models  are  not  correct  as 
approximations  to  the  optimal  (e.g.,  regulator)  control  laws  for  the 
original  problem.  While  control  and  filtering  theory  with 
homogenization  is  not  very  advanced  at  this  stage,  it  is  nevertheless 
possible  to  analyze  some  prototype  problems  to  a  point  where  the  basic 
features  of  the  theory  are  clear.  Our  work  has  contributed  to  this 
process,  but  much  more  needs  to  be  done. 


1.  Background:  dynamical  Control  of  Flexible  Structurea 

Interaction  of  the  control  syateo  with  the  structural  dynjuaics  of 
the  nechanical  system  is  one  of  the  fundamental  issues  in  large  space 
structure  applications.  Our  work  is  intended  to  contribute  directly 
to  understanding  this  interaction  by  using  models  which  capture  the 
essential  distributed  character  of  the  system,  and  using  analytical 
techniques  which  preserve  the  character  of  the  physical  system  in  the 
model  simplification  process.  The  methods  we  have  used 
¥iener-Hopf/spectral  factorization  methods  for  design  of  distributed 
control  systems  and  homogenization/ asymptotic  analysis  for  model 
simplification  --  have  tremendous  potential  for  the  smalytical 
treatment  of  complex  structural  control  problems,  including  the 
synthesis  of  computer'-aided-design  methods  for  large  space  structures. 
In  Phase  I  of  this  project,  we  have  concentrated  on  the  treatment  of  a 
few  simple  prototype  systems.  The  methods  may  be  adapted  to  treat 
complex  structures.  They  do  not  require  substantial  extensions  for 
such  cases.  Rather,  computational  algorithms  which  translate  their 
strengths  into  effective  design  tools  need  to  be  developed. 

The  main  emphasis  in  the  first  part  of  this  work  has  been  the 
adaptation  and  enhancement  of  certain  Viener-Hopf  methods  for  control 
system  design  used  by  J.  Davis  for  the  treatment  of  linear,  dynamic, 
distributed  parameter  models  of  flexible  structures  (Davis  1978, 
1979a,b  1982)  (Davis  and  Barry  1977)  (Davis  and  Dickenson  1983). 
Davis  and  his  colleagues  developed  a  frequency  domain  methodology  for 
computing  optimal  (regulator)  feedback  gains  for  linear  distributed 
parameter  control  systems  by  Viener-Hopf  spectral  factorization.  The 
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numerical  algorithms  for  executing  the  spectral  factorization  were 
based  on  some  earlier  work  of  P.  Stenger  (1972).  ¥e  have  adapted  the 
Davis-Stenger  methodology  to  the  problem  of  vibration  control  of 
flexible  structures.  A  generic  problem  of  this  type  is  the  figure 
control  of  a  large  apace  antenna.  Ve  have  carried  out  the  analysis 
and  computed  the  optimal  feedback  regulator  control  laws  for  several 
examples  including  the  Euler-beam  and  a  two-dimensional  prototype 
(experimental)  system  studied  by  J*  Lang  and  D.  Staelin  (Lang  and 
Staelin  1982a,b). 

7his  portion  of  the  research  has  demonstrated  the  effectiveness 
of  frequency  domain  —  spectral  factorization  methods  for  the  design 
of  control  and  state  estimation  algorithms  for  flexible  stiructures 
described  by  linear  distributed  parameter  models  (hyperbolic  partial 
differential  equations).  In  this  approach  it  is  not  necessary  to 
reduce  the  models  to  finite  dimensional  (lumped  parameter)  models  at 
the  outset  of  the  design  procedure.  The  infinite  dimensional 
character  of  the  system  is  preserved  throughout.  The  spectral 
factorization  methodology  avoids  the  difficult  numerical  problems 
associated  with  the  solution  of  the  Riccatl  partial  differential 
equations  which  arise  in  the  time  domain  approach  for  designing 
stabilizing  controllers.  In  this  way  distributed  phenomena,  like 
travelling  waves,  which  characterize  the  macroscopic  dynamics  of 
flexible  structures  are  retained  in  the  model,  and  their  interaction 
with  the  control  system  is  preserved  in  the  analytical  design  process. 
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1 . 1  Generic  Hodels  of  the  Qynaaice  of  Flexible  Structures 

The  flexible  structures  treated  in  this  work  are  assumed  to  be 
continue  described  generically  by  the  system  of  partial  differential 
equations 

(1)  m(x)h^^(t,x)  +  Dgh^Ct^x)  +  AQh(t,x)  -  F{t,x) 

where  h(t,x)  is  an  n-vector  of  instantaneous  displacements  away  from 
its  equilibrium  of  the  structure  S,  a  bounded  open  set  in  with 
smooth  boundary  S«  The  mass  density  m(x)  is  positive  and  bounded  on 
S.  The  damping  term  contain  both  (asymmetric)  gyroscopic  and 

(asymmetric)  structural  damping  effects*  The  internal  restoring  force 
term  A^h  is  generated  by  a  time- invariant  differential  operator  Aq 
specific  to  the  flexible  structure*  For  most  cases  of  interest,  Aq 
may  be  taken  to  be  an  unbounded  operator  with  domain  1>(Aq)  containing 
smooth  functions  with  the  appropriate  boundary  conditions  which  is 
dense  in  the  Hilbert  space  ■  I«^(S)  with  the  natural  inner  product, 
<*,*  >*  In  many  cases  Aq  has  a  discrete  spectrum  with  associated 
eigenfunctions  which  constitute  a  basis  for  L  (S)* 

The  applied  force  distribution  F(t,x)  generally  has  three 
components 

(2)  F(t,x)  -  F^(t,x)  +  F^(t,x)  ♦  Fj^(t,x) 

where  (t,x)  is  a  vector  of  exogeneous  disturbance  forces  and 
torques,  F  (t,x)  is  a  continuous,  distributed,  controlled  force  field 

w 

(as  in  an  electrostatically  controlled  system);  and  F2^(t,x) 


represents  the  control  forces  due  to  discrete  actuators 


(3) 


(t,x) 


.5^  bj(x)u^(t)  ^  BQU(t) 


The  actuator  amplitudes  are  u  (t)  and  the  actuator  influence  functions 
bj(x)  are  typically  elements  of  (which  usually,  but  not  always, 
approximate  delta  functions  6  (x-x^)*  Observations  are  usually 
assumed  to  arise  from  a  finite  number  p  of  sensors 


(4) 


yi(t) 


^*=3 'Vo 


j»i , . - 


or 


(4') 


y  (t) 


y(t) 


where  the  position  and  velocity  influence  functions  c^ ,  c^ , 
J  •  1,...,p  are  elements  of  H  which  may  represent  point  devices. 
Note  that  Bq:r”  H  ,  sH®  »  «iid  CjjtH  -*■  H  are  bounded. 


The  control  problem  for  (l)-(4)  is  to  choose  the  discrete  control 
amplitudes  u^Ct),  J  ■  l,...,m,  and  the  distributed  control  forces 
F_(t,x),  based  on  the  observations  y.(t),  i  >  l,...,p  to  maintain  the 

C  X 

state 


as  close  to  its  equilibrium  position  (nominally  zero)  as  possible.  If 
the  disturbances  are  transient,  this  nay  be  accomplished  by  using  a 
regulator  control  law  which  minimizes  the  quadratic  performance  index 


(6)  J  ta)  -  /q  [  q  tv.v)  +  au^  (t)u  (t)  +  Q  (P^,  r^)  ]dt 


where  q,  Q  are  non-negative  quadratic  forms,  and  a  is  a  positive 
parameter.  This  is  the  generic  control  problem  surveyed  in  (Balas 
1982).  It  includes  boundary  and  interior  control  of  vibrating 
strings,  moibranes,  thin  beams,  and  thin-plates. 


1.2  State  space  models  and  modal  control 

Suppose  for  the  moment  that  is  oymmetric  with  compact 
resolvent  and  spectrum  bounded  from  below.  The  spectrum  of 
consists,  therefore,  of  isolated  eigenvalues  , 


(7)  ^  Xj  ^  •  •  • 

and  eigenvectors  such  that  .  Assume  X^  >  0.  Then  A^ 

satisfies 


(8)  <AQh,h>Q  S  ellh||J,e>  0 

1  M 

and  Aq  has  a  square  root  A’  .  Let  D(A^)<^L^(s)  be  the  domain  of 
Aq,  and  D(A  J  )c:l2(s)  be  that  of  aJ  .  Let  H  -  L^S)  x  L^S),  and 
consider 


,  [vl  fw  "I  i 

'5)  A  «  ,  veD  (Aq)  ,  weD  (Aq) 

w  tvJ 


(10)  B 


C-[Cq  Cq] 


80  that  B:R”  x  h”  -►  H  and  C;H-^  R**.  With  v(t)  defined  by  (5)  we  have 


(11)  ^  V  (t)  -  Av  (t)  +  Bu  (t) ,  y  (t)  »  Cv  (t) ,  V  (O)el^ 


5 


whlcA  is  the  state  space  description  of  the  original  problem  with  the 
additional  assumptions  on  A. 


The  energy  inner  product  <*  ,•  >  defined  on  IL  is 


And  so,  in  the  energy  norm  ve  have 


(13)  +  <h^,h^> 

which  is  a  measure  of  the  total  potential  and  kinetic  energy  in 
(h,h^).  The  operator  A  on  (l^,  <*  generates  a  unitary  group  U(t) 

(Treves,  1975)  and 


(14)  0(t)VQ 


I 

k-1 


-1  . 

cos  uij^t  0)  sin  Uj^t 

*k‘°> 

sin  Uj^t  cos  uij^t 

for  any 


(15)  ^0  " 

00 

z 

k-1 

ak(0) 

♦k  '  "i 

Thus,  when 

u(t)  =  0 

in 

( 1 1 ) ,  energy  is 

!lu(t)v||  - 

1 1  „  1 1 
"^o"e' 

for  any  Vq  e  H  .  For  any 

differentiable 

,  the  solution  of  ( 1 1 )  is 

(16)  v(t) 

-  U(t)VQ  + 

/q  0(t-T)Bu(T)dT 

In  fact. 


(17)  v(t) 


£ 


with  [i^(t),  k  -  1,2,.*.,  dwfinad  by  ordixiary  differential 

equationa . 

)c 

Sy  introducing  finite  dimensional  subspaees  H  -  span 
k  >  1,2,...,k}  of  ,  one  can  construct  finite  dimensional  modal 
approximations  to  the  system  (11);  and  from  these,  finite  dimensional 
control  problems  whose  solutions  may  be  used  to  compute  suboptimal 
control  laws  for  the  infinite  dimensional  control  problem  defined  by 
(6)(11).  (See,  e.g..  Bales  1978,  1982).  The  feedback  controls 
obtained  in  this  way  will  stabilize  the  first  K  modes  of  the 
distributed  system.  However,  as  noted  in  (Bales  1978)  in  all  but  a 
few  special  systems,  the  ccntrol  actions  will  excite  the  higher  order 
modes.  This  "spillover"  effect  invariably  degrades  system 
performance.  While  this  phenomenon  has  received  considerable 
attention  in  the  literature,  it  is  the  unavoidable  companion  of  design 
methods  based  on  lumped  parameter  models. 


In  this  research  we  take  a  different  approach  to  the  control 
system  design  which  deals  directly  with  the  infinite  dimensional 
system.  The  method  uses  a  frequency  domain  formulation  of  the  control 
problem,  analogous  to  the  setup  for  finite  dimensional  problems  in 
(Willems  1971)*  and  a  spectral  factorization  algorithm  to  compute  the 
feedback  gain.  The  formal  algorithms  are  described  in  section  2. 
Before  developing  the  mathematics  it  is  useful  to  look  at  some 
examples  and  prototype  systems  which  illustrate  the  basic  features  of 
the  control  problem. 
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1 . 3  Control  of  a  vibrating  flexible  string 


Small  vibrations  of  a  flexible  string  may  be  described  by 
(16)  *  r(Oh.  - 

where  p(z)-  >  Pq>  0  is  the  linear  mass  density,  p(z)  >  Pq>  0  is  the 
modulus  of  elasticity,  and  we  assume  that  p  ,  p,  q,  r  are  twice 
continuously  differentiable.  Suppose  that  the  space  ze[0,l]  and  time 
t  have  been  normalized  to  dimesionless  coordinates.  The  system  can  be 
put  in  a  standard  form  by  changing  ihe  independent  variable 

X./  )*as 

with  L  =  x(l )  we  have 


®  h  +  b  (x)  h  =  0 

t  X 

0  ^  t 

The  coefficients  a(x)  and  b(x)  are  continuously  differentiable 
functions  of  x.  Defining 


(19) 


tt 


-  h  + 

XX 

0  ^  -  T. 


(20)  V  =  h^(t,x),  w  =  h^(t,x) 


we  have 


The  appropriate  boundary  conditions  are 

0£QV(t,0)  +  0QW(t,O)  =  0 

(22)  a^.v(t,L)  +  S^w(t,L)  =  u(t) 


8 


with  the 


conditions  J*  tl»  tl  inposed  to  avoid 

pathologies.  Here  6g  **  0  corresponds  to  a  fixed  endpoint,  while 
-  0  permits  an  end  to  move  freely  along  the  h  axis,  and  0, 

0  represents  an  end  free  to  move  but  with  positive  or  negative 
friction.  The  function  u(t)  is  a  boundary  control. 

One  can  generalize  (21 )  to 


and  the  real  coefficients  are  continuously  differentiable  on 
0  S  xS  L.  By  studying  the  finite  time  controllability  of  (22)-(24) 
D.L.  Russell  (1972)  was  able  to  prove  some  interesting  properties  of 
the  eigenvalues  and  eigenfunctions  of  A.  In  particular,  if  the 
(complex)  eigenvalues  of  A  are  then  {  e^J«,*k  -  1,2  ...}  forms  a 
Riesz  basis  for  the  space  L^O,  2L].  Moreover,  there  is  a  unique 
control  u  e  L^[0,T]  T*‘2L  (recall  that  all  variables  are  dimensionless) 
which  takes  the  solution  of  (22)-(24)  to  zero  at  t“T“2L  from  arbitrary 
initial  conditions  (in  the  space  L  ) 


(25)  v{0,x) 


v  (x)  ,  w(0,x) 


w  (x) ,  0  -  X 


< 


L 


(26)  k/J  |vq(x)|^  +  |wQ(x)|^dx  ^  /q  |u(t)|^ 

SK/L1^o(x)|2  +'  |uj,(x)|2  dx 

for  some  positive  constants  k,  K.  The  time  T  =  2L  is  "critical".  In 
general,  it  is  not  possible  to  make  the  transfer  in  T<2L;  and  for 
T  >  2L  there  vill,  in  general,  be  many  controls  which  accomplish  the 
transfer.  ^y  considering  the  special  spectral  structure  of  the 
operator  A  and  its  adjoint  h*,  Russell  was  able  to  show  that  the 
unique  control  u(t),  0  it  i2L,  accomplishing  the  transfer  could  be 
synthesized  by  a  bounded  linear  functional  of  the  state  in  (21). 

From  this  analysis  it  follows  that  the  optimal  regulator  problem 
for  (21),  that  is,  the  problem 


I 


k 


t 


V*. 


(27)  min  /.  tu  (t)  +  {/  (t,x)|  +  |w(t;x)l  dx)]dt 

ueu  . 

ad 

subject  to  (21)  (22)  with  admissible  controls  consisting  of  bounded 
(linear)  feedback  functionals  of  the  state  has  a  unique  solution  which 
produces  a  finite  optimal  cost.  The  problem  (25)-(25)  (27)  is  the 
simplest  example  of  the  class  of  control  problems  treated  in  this 
work.  It  is  a  one  dimensional  version  of  the  two  dimensional 
prototype  system  discussed  next. 

1.4  Control  of  a  two-dimensional  hyperbolic  system 

In  an  interesting  paper  J.  Lang  and  D.  Staelin  (1982b)  studied 
the  dynamical  control  properties  of  a  simple  experimental  system  as  a 
prototype  of  an  antenna  design  using  electroatatic  control  to  maintain 
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.-•'.'vV 


••  .v . 


.'•3 


t: 


I 


■ 


the  antenna  shape  (Lang  and  Staelin  1982a).  The  experimental  system 
consisted  of  a  flexible,  conducting  wire  mesh  (about  1  m  )  suspended 
vertically  in  tension  by  rigid  rectangular  boundaries  and  biased  by  a 
high  voltage  source.  A  parallel  surface  of  equal  dimension,  spaced  a 
short  distcmce  normal  to  the  mesh,  supported  a  3  x  3  array  of  fixed 
conducting  plates  independently  addressable  through  bipolar,  variable 
low  voltage  sources  which  collectively  served  as  a  distributed, 
electrostatic  control.  A  similar  set  of  plates,  equally  spaced  from 
the  mesh  on  the  opposite  side,  served  to  capacitlvely  sense  mesh 
deflections.  The  balancing  electrostatic  pressures  on  the  mesh 
produced  a  grounded-control  equilibrium  geometry  in  which  all  three 
surfaces  were  parallel. 

A  regulator  control  law  was  designed  to  modulate  the  voltages  on 
the  9  actuator  plates  in  response  to  (filtered)  observations  of  the 
mesh  deflections  from  the  9  sensor  plates.  Finite  dimensional  modal 
models  representing  the  dynamics  of  1  to  3  modes  were  used  in  the 
control  system  design.  The  basic  linear  -  quadratic  -  Gaussian 
regulator  control  law  was  not  satisfactory  in  certain  experimental 
regimes  (high  bias  voltage)  due  to  unmodeled  physical  factors. 
Modifications  were  necessary  to  achieve  mesh  stabilization  in  these 
cases.  Spillover  effects  were  also  observed  and  compensated. 

In  (Lang  and  Staelin  1982b)  the  mesh  was  modeled  as  a  flexible 
membrane  in  tension.  The  transverse  mesh  deflection  h(t,x,y),  defined 
as  positive  toward  the  sensor  plates,  satisfies 

«  *  V„  *  Vyy  -  Oh^  *  t 


(23)  Hh 


Here  M  is  (uniform)  mass  density,  llx  •  ^6  (uniform)  coefficients  of 
mesh  tension,  and  0  is  a  viscous  damping  coefficient.  Assuming  a  long 
wave  model  for  the  electric  field  between  the  mesh  and  plates,  the  net 
transverse  electrostatic  pressure,  f(t,y,7),  acting  on  the  mesh 
satisfies 


I  t 


f  =  i  e. 


(V-u)^l 
o  (H-h)  (H+h)  I, 


where  u  »  u(t,x,y)  is  the  potential  of  the  actuator  plates,  V  is  the 
mesh  bias  potential,  H  is  the  electromechanical  plate- to-mesh 


I  I 


i  E 


separation,  and  is  the  permittivity  of  free  space.  Assuming 
lhl«  H  and  iti|  «  V 

equation  (29)  was  linearized,  and  the  resulting  linear  control  system 
studied  in  (Lang  and  Staelin  1982b)  was 

(30)  MN:t  =  Vxx  +  +  Kh  -  Bu(t) 


where  K  -  2^V^/H^and  B  -CqV/H^.  Defining  S  -  [O,  ]  x  [o,  Lg]  as 
the  location  of  the  mesh,  the  mesh  boundary  conditions  required  zero 
deflection  at  the  perimeter. 


t  * 


r  V. 
> 


I  t 

I 

I 

I  c 


If  we  identify  Ah  as  the  linear  operator  on  the  right  in  (30), 
then  the  eigenfunctions  of  A  are 

(51)  ”  sin(mTTx/L^)  .  sin  (niTy/Lg) 

0  ^  X  ^ 


and  the  corresponding  eigenvalues  are 


•■-vN 

vS-’i 


.I> 


r 


*“a  P 


These  define  the  "open- loop"  natural  frequencies  of  the  systea. 
Notice  that  the  (a,  n)-node  is  unstable  if 


(35)  > 


2  e. 


Therefore,  if  V  is  large  enough,  a  finite  nuaber  of  aodes  are 
open- loop  unstable. 


The  experiaental  systea  in  (Lang  and  Staelin  1982b)  has  noise  in 
both  actuator  and  sensor  systeas.  This  noise  ras  represented  by 
Gaussian  white  noise.  The  overall  perforaance  index  used  in  the 
design  was 


(34) 


lia  e{— — — 

I 


,  (k+l)  T.  ^^B 

''kT  ^0  ■'o 


q  h  (t,x,y) 


+  (t,x,y)  dx  dy]dt} 


where  r  -  ,  with  Vq  the  dyncunic  range  on  the  voltage  control 

systea,  and  q  "  2  V/HVq  . 

The  study  of  this  saall  systea  provides  a  great  deal  of  useful 
inforaation  on  control  probleas  that  can  be  expected  for  certain 
classes  of  flexible  structures.  The  control  systea  perforaance 
reported  by  Lang  and  Staelin  provides  one  of  the  few  available 
benchaarks  against  which  alternative  control  systea  designs  aay  be 
tested.  ¥e  shall  consider  this  systea  further  in  section  4. 
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1 . 5  Control  of  a  simply  supported  beam 


The  Euler-Bemoulll  equation  for  the  dynamics  of  a  simply 
supported  beam  Is 


Mh 


tt 


Dh 


txx 


Elh 


(55) 


0  -  X  -  L 
h(t,0)  - 


xxxx 
0  -  t 


f (t,x) 


h(t,L) 


h  (t,0) 

XX 


h  (t,L) 

XX 


where  M  Is  the  mass  density  (per  unit  length),  S  Is  a  damping  ratio,  E 
Is  the  modulus  of  elasticity,  I  Is  a  moment  of  Inertia,  and  f  Is  an 
applied  force  distribution. 


If  we  Ignore  the  damping,  0-0,  and  normalise  other  parameters 

to  unity,  then  the  mode  shapes  -  eigenfunctions  are  ^j((x)  **  slnkicx  and 

2 

the  eigenvalues  are  (icrr)  .  Control  of  vibrations  of  the  beam,  may 
be  based  on  the  performance  Index 

(36)  j{u)  -  -Tq  *  qJh^Ct.x)  +  q^h2(t,x)dx)]  dt 


Numerical  studies  of  this  problem  were  reported  in  (Balas  1978).  A 
point  actuator  and  a  point  sensor 

(37)  f  (t,x)  =  u(t)  5(x  -  x) 
y  (t)  «  h  (t,x)  5  -X  ) 

were  used  to  effect  control  in  the  problem.  Spillover  into  the 
uncontrolled  residual  modes  produced  instability  In  the  simulations. 
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2.  Wiener-Hopf  Methods  for  Control  System  Design: 
The  Davis-Stenger  Algorithm 


The  connections  among  least  squares  optimization,  spectral 
factorization  and  algebraic  Riccati  equations  have  been  considered 
important  in  control  theory  for  many  years.  (See,  e.g.,  Anderson 
(1967),  Brockett  (1970),  Willems  (1971),  Molinari  (l973a,b),  Helton 
(1976),  and  the  references  therein.)  To  see  how  the  connection  arises, 
consider  the  standard  finite  dimensional,  infinite  time  linear 
regulator  problem 


r 

min  /  „ 
u  0 

+  |y(t) 

l^dt 

x(t)  = 

Ax  (t) 

+  Bu(t) 

,  X 

y(t)  = 

Cx(t) 

,  t  ^  0 

Suppose  A  is  a  stable  matrix  and  (A,B,C)  is  a  minimal  finite 
dimensional  triple  realizing  the  transfer  function 
(2)  G(s)  =  C(ls  - 


Then  the  optimal  feedback  control  for  (l)  is  given  by 
(3)  u(t)  =  -B*Kx(t) 

with  K  the  unique  positive  definite  symmetric  solution  to  the 
algebraic  Riccati  equation 


(4) 


A*K  +  K(A  -  BB*K) 


-C*C 


An  integral  equation  for  the  optimal  feedback  gain  may  be  derived 


from  (4)«  Let  s  be  a  complex  number  not  in  the  spectrum  of  -A*, 
a(-A*),  nor  in  the  spectrum  of  A-BB*K,  then 

K[  Is  -  (A-BB*K)r^  +  (-Is  -  A*)“^K  =  (-Is  -  A*)"^C*cC-IS  -  (A-BB*K)] 

t  -  \ 

From  stcmdard  results  (Brockett  1970), a  (A-BB*K)  is  contained  in  the 
open  left  half  of  the  complex  plane;  and,  by  assumption, a  (-A*)  is  in 
the  open  right  half  plane.  Let  F  be  a  closed  rectifiable  contour 
encircling  a (A-BB*K)  in  the  positive  sense,  and  integrate  (5)  along T . 

Since 

^  -1 

/pi  Is  -  (A  -  BB*K))  ds  =  I 

21T  i 

(6) 

— -  /r  “  A*  j”^ds  «*  0 

2Tf  i  '■ 


we  obtain 

1  -1  -1 

(7)  K  -  -  /pl-ls  -  A*1  C*C  I  Is- (A-BB*K)  J  ds 

2lTi 

and  so 

(8)  KB  -  -  /p  [  Is  -  A*]  "^0*0  [Is-(A-BB*K)]"^Bd3 

2iri 

Since  the  integrand  is  the  product  of  two  rational  functions,  the 
contour  may  be  deformed  to  yield 
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VjcVjkJT.*.  O-*.. 


:v 

C 


r>  ■■■ 


KB  »  -  ^UliW- A*]“^C*C  tiiw-  (A-BB*K)l"^Bda) 


The  spectral  factorization  identity  (Brockett  1970,  Villens  1971) 


Httw)  -  I  +  G*  ttw)  G  aw) 

=  I  +  B*  ( -liw  -  A*J  ~^C*C  [  liw  -  A] 

*  F*  aw)  F  aw) 

-  [I  +  B*  (-liw  -  A*)‘^  KBl  [I  +  B*  K  ffiw  -  A)  S] 


and  the  identity 


ai)  C  [liw  -  (A  -  BB*K)]"'S  -  Cttiw  -  A)"^b[  I  +  B*K  ttiw  -  A)"S] 


when  used  in  (9)  gives  the  result 


(12)  ®**‘  “  —  ^lF*(iW)  1  '■H*R*(iW,A)C*CR(iw,A)dw 

2  IT 

where  R(s,A)  ■  [ls-A]”^is  the  resolvent  of  A. 


Hence,  to  compute  the  optimal  feedback  gain,  we  can  either  solve 
the  nonlinear  algebraic  Rlceati  equation  (4)  for  its  unique  positive 
definite  solution  or  we  can  carry  out  the  spectral  factorization  of 
I>G'*G  in  (10)  and  then  compute  the  integral  (12).  In  finite 
dimensional  control  problems  there  may  be  little  reason  to  favor 
formulation  of  the  computational  problem  in  one  setting  -  the  Riccati 
equation  -  over  the  other  -  spectral  factorization.  In  infinite 
dimensional  problems,  however,  the  spectral  factorization  method 
appears  to  have  superior  numerical  stability  properties  over  direct 
integration  of  the  Riccati  equation.* 


•  e 


'  s  o. V'. V 
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2«1  Davis's  Method 


In  a  series  of  papers,  J.  Davis  and  his  students  (Davis  and 
Barry  (1977),  Davis  ( 1978a, b)  (1979)  (1982),  Davis  and  Dickinson 
(1983))  have  explored  the  application  of  spectral  factorization 
methods  for  control  system  design  to  a  class  of  distributed  parameter 
models  of  long  trains  with  multiple  locomotives.  The  control  problem 
is  to  modulate  the  acceleration  of  individual  locomotives  to  minimize 
deviations  in  coupler  stress  throughout  the  train.  Disturbances 
Include  passage  of  the  train  over  a  grade,  which  tends  to  set  up  a 
"travelling  wave"  along  the  train  of  stress  deviations  from  nominal. 
The  first  approach  to  this  problem  which  comes  to  mind  is  to  write  out 
the  equations  of  motion  of  the  cars  and  locomotives  in  the  train  and 
formulate  an  optimal  control  problem  for  the  overall  system.  The 
large  dimension  of  the  resulting  model  and  the  absence  of  cmy  special 
structure  Inhibits  this  approach.  Decentralized  control  schemes 
(HcLane,  Peppard,  and  Sundareswaran  (1973),  Gruber  and  Bayoumi  (1982)) 
are  not  particularly  effective  for  these  problems.  As  Davis  and  Barry 
(1977)  have  observed,  aside  from  the  difficulties  in  solving  large 
scale  control  problems,  one  has  trouble  estimating  the  effects  of 
system  parameter  changes  or  of  variations  in  the  number  of  units  in  a 
block  based  on  lumped  parameter  models. 

Davis  recognized  that  the  mass-spring  nature  of  the 
Interconnected  system  could  lead  to  traveling  wave  phenomena  setup  by 


"competing"  local  controllers  (locomotive  accelerations).  He  reasoned 
that  the  macrosopic,  widely  coupled  motions  of  the  elements  would 
contain  the  bulk  of  the  energy  of  motion  of  the  total  system.  This 
hypothesis  suggests  that  a  control  scheme  designed  to  suppress  such 
motions  would  achieve  substantial  reductions  in  the  coupler  stress 
levels. 


To  represent  the  system  in  a  fashion  which  would  capture  wave 

phenomena  most  naturally,  Davis  reformulated  the  system  as  a 

distributed  parameter  system  with  boundary  controls.  (Davis  and  Barry 

1977).  The  resulting  model  proved  to  be  mathematically  tractable. 

The  effects  of  changes  in  both  system  parameters  and  the  number  of 

units  in  a  block  were  readily  apparent.  In  fact,  an  increase  in  the 

number  of  elements  in  a  block  increases  the  validity  and  usefulness  of 

the  model.  In  contrast,  finite  dimensional  models  tend  to  become 
« 

increasingly  intractable  as  the  number  of  units  in  the  system 
increases . 


Davis’  modeling  technique  is  simple  and  instructive.  Consider 
the  mass  -  spring  -  damper  system  in  Figure  2.1.  The  dynamics  are 


"  S  *1  ■  ‘‘i-l  *  -  -i.!'  -  =  if  '”1-1  * 


where  x^(t)  is  the  deviation  of  the  1^  mass  from  its  nominal 
position.  The  continuum  approximation  to  this  system  may  be  developed 
as  follows:  Let  z  e  Co,l]  be  spatial  position  along  the  "rest 
length",  unity,  of  the  overall  system,  and  let  u(z, t)  be  the  deviation 
of  the  mass  at  rest  position  z  and  time  t.  Making  the  identification 


Figure  2.1  Discrete  element  model  for  viscoelastic  bar 


(14)  ua/N,t)  -v  X^(t),  i  -  0,1..... 


and  usin^;  the  numerical  differentiation  formula 


(^5)  2  [  u  (z+  2u  (z)  +  M  (z- 


1]  ^  a  u 


One  arrives  at  the  distributed  model  equation 

_3'*_v-#1\2  8u._,1,2  3u  • 

3t^  3z  3z  Jit 

as  an  approximation  for  the  motion  of  the  system.  Rewritten  in  the 

form  ^  3^u  3  , 

P  — r  *  —  s  (s.t) 

3t  3z 
^  K  If  . 


the  equation  describes  the  longitudinal  motion  of  a  visco-elastic  bar; 
the  term  s(z,t)  represents  the  stress  in  the  bar,  here  proportional  to 


a  linear  combination  of  the  strain  and  stress  rate  (Davis  and  Barry 
1977)  (Greenberg,  NacCamy  Mlsel  and  1968).  The  natural  boundary 


conditions  for  (17)  are  in  terms  of  8(z,t)  at  z*0  and  1,  i.e., 

(18)  s(0,t)  =  fgCt),  3(1, t)  -  fj^(.t) 

the  applied  forces  and  u(z,t),  u^(z,t)  at  t~0. 

Rescaling  time  at  t*  ■  t/(k/c),  the  system  (16)  may  be  rewritten 
in  dimensionless  form  as 


3  u  3  3  .  ,  3u  3  u  .  ^<*<1 

^  +  (_  +  ^)  o<x<l. 


3t‘ 


with  a  ■  c^/kHN^.  And  this  may  be  written 
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in  matrix  form  as 


0^x61,  O^t 


(20b) 


Vz  *  "ztl  "  ^0 
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(z,0),  u(z,0)  given 


Let  A  be  the  matrix  differential  operator  on  the  right  in  (20a). 
Defining  H  ■  L^[0,l]  x  L’[0,l]  with  the  energy  inner  product, 


(21 )  <  l"l  ,  I  S  1  >  *  J’o  ta'iP  +  vq)  dz 


then  A  on  H  has  domain 


D(A)  ■  ^ *  '*»v,v  absolutely  continuous  u  ,v  eL^  (0,1], 

(22)  ^  *  X  J« 

u  +  -  0  at  X  -  0,1} 

dense  in  H,  A  is  a  dissipative  operator,  and  A  is  the  infinitesimal 
generator  of  a  class  C  ^  contraction  s«Bigroup  on  H  (Davis  and  Barry 
1978,  Theorem  l).  Moreover,  the  resolvent 

(23)  R(s,A)  -  /  “  e"*'^T^  dt 

of  A  may  be  computed  explicitly*  (See  section  3  or  Davis  and  Barry 
1978.) 

It  is  not  possible  to  write  the  solution  of  (2)  in  the  strong 


(24)  -St 


where 


AU  +  Bf (t) 


0=  [  u^,  u  ]*  ,  f 
z  t 

H 


since  the  boundary  forces  correspond  to  generalized  function  "inputs". 
Using  L'^as  the  inverse  Laplace  transform,  Davis  and  Barry  treat  (20) 
in  the  form 


(26)  0(t)  -  U(0)  +  l“^  {G(s;z)  f  (z)} 
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vhere  6(8;z)  represents  the  "transfer  function  oatrlx"  associated  with 
the  boundary  value  problem  (20)* 

The  (optimal)  control  problem  involves  minimizing  variations  in 
the  stress  distribution  throughout  the  system.  The  quantity 


^  (*»t)  dz*  /q  Qi^(u^  +  dz  oi-a/n 


corresponds  to  the  total  stress  in  the  system.  Recalling  the  original 
approximation (13)- (16),  ve  have  the  correspondence 

(28)  /J  =82  ■>.  H  J  ^  ^  (2i  -  «i.i)  +  5  ,2^  -  )  1  = 

and  so,  the  natural  quadratic  cost  functional  is 

(29)  J-  ^  1^  +  1^  +  i  /J  (u^  +  u^^) ^  dz}dt 

This  formulation  includes  stress  contributions  from  spatial  modes 
of  all  wavelengths.  In  most  physical  systems  high  order  modes  will 
have  a  negligible  contribution  to  the  overall  behavior.  Using  to 

Xr 

denote  projection  onto  the  subspace  of  H  spanned  by  the  first  p 
eigenfunctions  of  A,  and  defining  the  system  output  as 

(30)  y  (2»t)  -  a[l,3^1  [TTpO]  U,t) 

the  final  formulation  used  by  Da'”’ls  and  Barry  (1978)  is 

7  {  |f„  (t)|=  ♦  |fj(t)|=  .i/lyfe.til^ax)  at 
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0(z,t)  -  0(0) (z)  +  {G(a»z)  f  (z)}  (t) 

y{z,t)  ■  all,3  ] [WO]  {z,t) 
z  p 

0(0)  e  D(A)C.L^[0,1]  X  L^[0,1] 

The  resulting  optimization  problem  is  a  distributed  control  problem 
with  state  cost  restricted  to  a  finite  dimensional  subspace. 

Davis  and  Barry  (1978)  compute  the  optimal  control  law  for  this 
problem  using  the  spectral  factorization  algorithm  described  earlier. 
A  key  step  in  the  procedure  is  application  of  a  numerical  algorithm 
for  spectral  factorization  due  to  F.  Stenger.  In  the  next  subsection 
we  summarize  Stenger 's  algorithm. 

2.2  Stenger 's  Algorithm  for  Spectral  Factorization 

To  evaluate  the  control  lav  for  a  given  problem  modeled  as  in  the 
last  subsection,  we  must  compute  the  spectral  factor  F(s)  appearing  in 
equation  (10).  That  is, 

(32)  F*(s)  F(s)  =  H(s)  -  I  +  G*(s)  G(s) 

where  G(s)  is  the  transfer  function  of  the  system  being  controlled. 
The  first  problem  is  to  determine  conditions  under  which  the  spectral 
factor  exists.  Since  G(s)  is  the  transform  of  a  real  vector  valued 
function,  which  we  assume  to  be  integrable  and  square  integrable,  it 
follows  that  S(a)  ■  G*(s)G(s)  is  a  Hermitian  positive  semidefinlte 

I 

(matrix  valued)  function  and  G*G  is  the  transform  of  a  function  which 
is  in  L^/1 
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convolution  algebra  I  •  L  ,  or  on  L  ,  by 

(37)  **+  L  f  (t)  e  -1+7  f  (t)  e 

0 

Stenger's  algorithm  is  used  to  provide  a  numerical  approximation  to 

the  causal  projection  operator.  Before  discussing  that,  we  note  that 

under  the  assumptions  on  G(s),  and  therefore  on  S(s),  that  the 

iteration  may  be  shown  to  converge  from  a  suitable  initial  guess  to 

2 

the  uniquely  defined  (in  L  )  causal  spectral  factor  F(s). 


The  algorithm  (36)  has  a  particularly  simple  form.  As  noted  in 
(Davis  and  Dickenson  1983),  the  computation  of  P^[*]  is  the  most 
difficult  step,  but  Stenger's  result  takes  care  of  that.  The 

numerical  approximation  in  (Stenger  1972)  takes  the  form 

» 

(38)  p*  Ifl  to)  -  f.  to)'v  J  Jod.  .  Jh)  ^ 

m 


Here  h  is  the  step  size  and  and  are  parameters  defined  by 

Stenger.  Specifically, 


(39)  a 


TTh 


1-iq 
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m 


4k  iK 


m  2 
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where  q  is  a  parameter  chosen  in  the  algorithm  and 


K-  (a/b)‘ 


k  -  iiTb 


The  step  size  Is  chosen  based  on  the  bandwidth  of  the  transfer 
function  f(s)«  If,  in  (38),  one  chooses  to  sample  the  projection  f^ 
at  the  same  sample  points  as  f,  then,  as  observed  in  (Davis  and 
Dickenson  1983) 

(41)  f.  (kh  +  ih)  '''  ^  f(jh  +  Jh)  E  _ ^ _ 

+  j  m  ■ 

(k-j)h-a  h  +  ih 
Si 

and  it  is  clear  that  the  required  calculation  is  a  convolution.  Since 
the  range  of  sample  points  is  finite  this  is  naturally  implemented  by 
a  fast  Fourier  Transform  (FFT);  and  this  substantially  improves  the 
computational  time. 

Since  we  must  compute  (F*)  S(F  )"^  in  the  iteration  (36)  for  the 

n  n 

spectral  factor,  it  is  beat  to  rewrite  the  iteration  as 

(42) 

and  execute  it  in  this  form.  As  observed  in  (Davis  and  Dickenson 

1983),  the  last  factor  in  this  expression  is  a  perturbation  of  the 

identity*  (since  (F*)“^S(F  )"^  -  I  O),  and  this  has  natural 

n  n 

advantages  in  the  numerical  realization  of  the  iteration. 

As  suggested  in  (Davis  and  Dickenson  1983),  a  suitable  choice  for 

the  Initial  guess  for  the  spectral  factor  F^  is  the  diagonal  matrix 

whose  elements  are  (scalar)  spectral  factors  of  the  corresponding 

diagonal  elements  of  S.  This  choice  implies  that  the  matrix 

(p*)“is(F  )“^  is  a  matrix  with  ones  on  the  diagonal  and  with  all  the 
n  n 

off-diagonal  elements  less  than  one  in  magnitude.  This  tends  to 
prevent  the  iteration  from  blowing  up.  The  diagonal  factors  may  be 


obtained  by  an  FFT  implementation  of  the  scalar  algorithm  in  (Stenger 
1972). 

The  method  as  described  here  vas  implemented  directly  on  the 
problem  of  controlling  the  dynamics  of  the  Buler-Bemoulli  beam.  The 
results  are  shown  in  the  next  section.  A  careful  consideration  of 
Stenger 's  method  suggests  an  alternative  implementation  of  the 
algorithm  which  takes  advantage  of  the  occurence  of  Hilbert  transforms 
in  the  course  of  the  computations  and  the  effective  use  of  these 
transforms  in  the  representation  of  the  causal  spectral  factor.  This 
observation  permits  an  efficient  numerical  realization  of  the  spectral 
factorization  procedure.  We  shall  develop  this  in  the  context  of 
design  of  stabilizing  controls  for  a  two  dimensional  flexible 
structure.  This  result  and  the  associated  algorithm  are  reported  in 


3>  Control  of  a  One-Dimensional  Structure 

The  control  of  simple  one  dimensional  structures  serves  to 
illustrate  the  general  analytical  methods  in  the  simplest  form.  One 
dimensional  models  can  also  represent  certain  components,  e.g., 
flexible  beams,  which  appear  in  composite  large  space  structures;  and 
they  may  represent  certain  two  or  three  dimensional  structures  with 
natural  symmetry.  In  this  section  we  consider  a  simple  system,  the 
Euler-Bemoulli  beam  in  detail,,  working  through  the  computation  and 
simulation  of  the  optimal  stabilizing  feedback  gain. 

3.1  Control  of  a  Flexible  Beam 

The  dynamics  of  an  undamped  flexible  beam  undergoing  transverse 
notion  are  described  by  the  Euler-Bemoulli  partial  differential 
equation 

m  u^^(t,x)  +  Elu^^  (t,x)  =  f(t,x) 

0<x<L,  t>0 

where  u(t,x)  is  the  transverse  displacement  of  the  beam,  f(t,x)  is  an 
applied  force  distribution,  m  is  the  mass  per  unit  length,  I  the 
moment  of  inertia,  E  the  modulus  of  elasticity,  and  L  is  the  beam 
length.  To  facilitate  comparison  of  our  results  with  earlier  work 
(Balas  1978b},  we  shall  assume  that  m,  E,  I,  and  L  are  all  unity.  The 
boundary  conditions  for  pinned  support  are 

u(t,0)  =  0  »  u(t,l) 

(2) 

u^(t,(a)  =  0  -  u^(t,l) 

The  beam  is  controlled  by  a  single  point  actuator 


(3)  f(t,x)  » -^  6(x-a)v(t)  ,  0<a<l 

and  a  single  sensor  measures  displacement 

(4)  y(t)  =u(t,b),  0<b<l 

The  system  is  deterministic  and  actuator  and  sensor  dynamics  are  not 
modeled. 


Balas  (1978b)  designed  a  feedback  controller  for  the  first  three 
modes  of  the  beam  which  minimized  the  (unweighted)  energy  in  those 
modes.  His  controller  includes  a  six-dimensional  Luenberger  observer 
to  reconstruct  the  state.  -  The  energy  in  the  fourth  (residual)  mode 
increases  rapidly  due  to  spillover  effects. 


Our  approach  to  this  problem  is  based  on  the  infinite  dimensional 
model.  Taking  the  Laplace  transform  of  (l)>  we  have 


'^xxxx  ®  U(s,x)  -  F(s,x) 


(5)  U(s,0)  =  0  -  0(s,l) 


U  (s,0) 

XXX 


0  (s,l) 

XX 


the  Green’s  function  for  (l)  (5)  satisfies 

(6)  Gj^(s,x|x^)  +  s^G(s,x|x'  )  -  5(x-x') 


Consider  G  in  the  form 


(7)  G(s,x|x'‘) 


[A  sinhv'S^x)  0<x<x' 

B  sinh  JT  (1-x) )  x'<x<l 


with  A  and  B  complex  constants  to  be  determined.  Note  that  the 


boundary  conditions  in  (5)  are  satisfied  by  (7).  At  x*x'  we  have 


X*  +  e 


x'  +  e 


G  dx  +  s/  Gdx  =  l 

,  xxxx  , 

x'  -  e  x'  -  c 


(See  Tai  (1971 )  •  From  (6)  (8  )  we  have 


x'  +  e 

(9)  G  I  =1 

XXX  '  , 

X*  e 


and  from  this 


A\  -1 


B  j  sinh  s  1  sinh  s  x* 


sinh  s  d-x') 


which  gives  the  Green's  function. 


The  transfer  function  relating  the  input  f(t,x)  in  (?)  to  the 
output  y(t)  in  (4)  may  be  written  down  immediately  from  6(s,z|x'). 


(11)  Y(s)  =  U(s,b)  =  /J  C{s,x|x')  F(s,x')dx'  =  -  G(s,b|a)  V  (s) 


with  V(s)  the  Laplace  transform  of  v(t).  Identifying 


(12)  (s)  =  G(s,bla) 


we  have 


(13)  (s)' 


-  (sinh  a  a)  sinh  (1-b)  s) 


s  8inh/s~ 


,  x'>a<b^ 


-  (sinh  Cl-a)  a)  feinh  b  s) 


yi” 8^' ^  sinh/T 


,  x“ei>b^ 


Balas  (1978)  uses  a  ,  b  ;  and  for  this  case  we  have 

6  6 


(t4)  Tgb)  .  T  (a)  - 

•g  T  2  a  sinh/ s 


To  use  the  Davis-Stenger  algorithm  as  described  in  section  2,  ve 
must  compute  the  spectral  factorization  of 

(15)  H(s)  =  1  +  T*  (s)  T_  (s)  =  F*  (s)F(8) 

B  B 

Substituting,  ve  have 

(16)  Httoi)  -  1  +  i  feinh /Ia)/6)^  (sinh /-iU)/6)  ^ _ 

(iu)  3/2  (-iu)  3/2  sinh  /iu)  sinh /  -iu 

and  the  computational  problem  reduces  to  computing  the  spectral  factor 
F* (s)  from  (^g ) ,  and  then  computing  the  optimal  stabilizing  feedback 


(17)  B  K  a  ijjp  /^[F  (iu)  ]  G  (iu))  CRCLu,  A)  du) 


by  numerical  Integration. 


3«2  Numerical  Results 


In  the  computations  and  simulations  which  follow,  we  consider  the 
same  hecun  parameters  used  by  Balas  (unit  length,  with  parameters 
normalized  to  one,  and  zero  damping).  The  numerical  requirements  of 
the  algorithm  are  not  changed  for  more  realistic  choices  of 
parameters.  In  addition  to  the  case  considered  by  Balas  (with  one 
controller  and  one  observer  at  opposite  ends  of  the  beam) ,  we  also 
consider  the  effect  of  increasing  the  number  of  controllers  and 
observers,  and  finally  the  effect  of  delays  in  the  control  loops.* 

Practical  Implementation  of  the  algorithm  requires  frequency 
sampling  of  the  transfer  function  and  spectral  factors,  and  spatial 
and  frequency  sampling  of  the  resolvent.  It  was  determined 
experimentally  that  for  the  given  beam  parameters,  a  frequency  range 
of  plus/minus  30HZ  is  adequate,  since  the  gain  from  input  to  output  in 
the  range  beyond  this  is  insignificant,  see  Figure  3.1.  The  figure 
also  indicates  a  very  smooth  dependence  of  the  transfer  function  on 
frequency,  so  that  the  256  sample  points  used  in  the  program  are  quite 
adequate.  Spatial  discretization  is  done  using  100  equidistant 
points. 

An  example  of  the  feedback  gain  is  given  in  Figure  3*2,  for  the 
velocity  state  variable,  and  in  Figure  3.3  for  the  displacement  state 
variable.  (See  Appendix  A  for  other  gains  corresponding  to  different 

•The  issue  of  the  effects  of  the  delay  on  the  control  action  and 
the  system  stability  was  raised  by  Dr.  J.  Burns,  formerly  with 
AFOSR,  and  now  with  VFIdSU,  Blacksburg,  VA.  We  are  grateful  to 
Professor  Bums  for  his  inputs  on  this  problem. 
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arraxigements  of  the  controllers  and  outputs.)  A  sharp  peak  at  the 
point  of  observation  for  the  deflection  state  indicates  the  control 
effort  to  reduce  the  deviation  at  this  point  to  zero,  since  only  this 
point  contributes  to  the  optimization  cost.  In  this  case  ve  have 
penalized  the  state  deviation  at  the  observation  points  (in  the  cost 
criterion)  500  times  more  than  the  control,  so  this  is  "cheap 
control." 

The  feedback  gains,  one  for  the  speed  state  and  the  other  for  the 
deflection  state  are  integral  operators  as  defined  by  (17).  These 
gains  are  computed  off-line.  Computation  of  the  input  function  for  a 
given  time  requires  evaluation  of  the  integral  operator  B  K  acting  on 
the  state.  This  is  accomplished  in  the  program  by  approximating  the 
Integral  as  a  sum  of  piecewise  constant  functions. 

In  controlling  a  physical  beam  one  would  need  an  observer  for  the 
deflection  and  velocity  variables  that  would  use  (point)  observations 
of  the  beam  deflection  as  Inputs.  For  the  simulation  results  here  the 
deflection  and  its  derivative  are  obtained  by  numerical  integration. 


V 


The  case  studied  by  Balas,  where  the  controller  and  the  observer 
are  at  the  opposite  ends  of  the  beam,  exhibits  poor  observability  and 
controllability,  which  is  reflected  in  the  large  control  efforts  and 
long  settling  tines  to  stabilize  the  notions  (see  Figure  3.6).  Using 


il 


C 


more  observers  and  controllers  substantially  improves  the 
"controllability"  of  the  system.  For  example,  in  the  case  of  three 
collocated,  equidistant  controllers/measurements,  the  margin  of 
improvement  can  be  seen  by  comparing  Figure  3.5  with  Figure  3.6.  See 
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also  Appendix  A  for  the  time  responses  corresponding  to  these  two 
cases* 

Belays  in  physical  control  loops  are  inevitable  due  to  the  finite 
time  necessary  to  process  measurements  and  compute  the  resulting 
controls.  Analytical  treatment  of  delays  using  time  domain  models  is 
not  nearly  as  convenient  as  it  is  in  the  frequency  method  described 
here.  For  example,  if  the  delay  is  T,  than  we  need  only  multiply  each 
element  of  matrix  H(s)  in  (15)  by  a  factor  exp(-sT).  Therefore,  HCs) 
is  invariant  with  respect  to  the  delay,  and  the  critical  part  of  the 
gain  computation  algorithm,  i.e.,  spectral  factorization,  need  not  be 
recomputed  for  the  delayed  case.  Of  course,  the  delay  appears  in  the 
transfer  function,  and  so,  a  new  gain  must  be  computed  for  each 
different  delay.  (Numerical  results  displaying  effects  of  delays  are 
discussed  later.) 

To  validate  the  program,  we  simulated  the  response  of  a  beam 
subject  to  an  initial  disturbance  in  the  form  of  one  of  the  spatial 
modes,  and  with  a  feedback  based  on  the  optimal  control  discussed 
above.  Figure  3.5  is  an  example  of  a  time  variation  of  the  beam  at 
the  observation  point.  (Other  such  examples  are  given  in  Appendix  A.) 
Note  that  without  the  control  this  response  would  represent  an 
undamped  oscillatory  motion  since  the  beam  model  does  not  include 
damping.  Damping  has  a  stabilizing  effect  in  this  system;  and  the 
control  action  is  enhanced  if  damping  is  present. 


The  plots  suomarising  nufflerlcal  work  on  the  example  problem 
appear  in  several  groups  in  Appendix  k,  and  they  are  grouped  as 
follows.  For  a  given  delay,  an4  a  given  number  and  position  of 
controls  and  observations,  we  first  display  gains  for  deflection  and 
velocity  states,  both  for  each  of  the  inputs.  Hext  is  a  group  of 
plots  showing  the  beam  response  to  the  optimal  control  when  the  beam 
is  initially  displaced  in  the  form  of  one  of  the  first  three  spatial 
modes.  (Significant  components  of  the  matrix  H(i  )  are  well  below 
30Hz,  i.e.,  below  the  frequencies  induced  by  the  4th  spatial  mode.) 
The  meaning  of  the  plot  titles  is  as  follows:  "Beam  at  x,  yth  mode,” 
means  that  the  beam  is  initially  displaced  by  the  y-th  spatial  mode, 
and  the  deflection  of  the  beam  is  observed  as  a  function  of  time  at 
point  X.  On  the  control  plot  we  Indicate  the  position  of  the  point 
control  and  the  time  evolution  of  the  control  at  that  point. 


While  the  main  purpose  of  conducting  these  experiments  was  to 
verify  the  control  algorithm,  several  phenomena  were  observed  from 
these  experiments.  Comparing  responses  of  higher  order  modes  with 
those  of  the  lover  order  modes,  it  is  evident  that  more  energy  is 
needed  to  control  higher  order  modes  (note  that  our  model  has  zero 
damping).  This  reflects  the  poor  controllability  and  observability  of 
the  higher  order  modes.  Second,  the  gains  for  the  deflection  state 
have  pronounced  peaks  at  the  observation  points,  suggesting  use  of 
localized  -  decentralized  feedback.  Unfortunately,  the  speed  gains 
indicate  much  more  spatial  '  coupling,  suggesting  that  decentralized 
control  schemes  may  not  be  effective  (at  least  for  the  parameter 
ranges  used  in  this  problem).  Third,  the  delay  has  a  substantial 


effect  on  the  performance  of  the  control.  (Compare  any  case  with 
delay  from  Appendix  A  with  a  case  with  no  delay.)  Nevertheless,  the 
stability  margin  is  remarkably  wide.  An  example  from  Appendix  A 
indicates  that  a  delay  equal  to  one  half  period  of  the  highest  mode  in 
the  chosen  spectrum  does  not  destabilize  the  system. 

The  numerical  results  presented  above  affirm  the  Davis-Stenger 
algorithm  as  a  practical  tool  for  vibration  control  of  flexible 
structures,  represented  here  by  an  Euler  beam  model.  The  results  of 
this  algorithm  provide  a  means  for  assessing  effects  of 
controller/observer  placement  on  the  system  performance,  as  well  as 
give  stabilizing  feedback  gains,  once  the  controller  locations  have 
been  selected.  It  was  also  demonstrated  that  the  underlying  model 
allows  an  efficient  treatment  of  delays  in  the  control  loop. 
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Figure  3-5-  Time  response  of  tbe  beam;  three  actuators 


4«  Control  of  a  Two-Dimensional  Prototype  System 


In  this  section  we  adapt  our  frequency  domain  control  system 
design  procedure  to  treat  a  prototype  two  dimensional  flexible  system 
-  a  membrane/mesh  whose  dynamical  behavior  is  sensed  and  controlled  by 
electrostatic  forces.  The  model  la  patterned  after  an  experimental 
system  studied  by  Lang  and  Staelin  (l982b}  as  a  paradigm  for  an 
electrostatically,  controlled  large  aperature  reflecting  satellite 
antenna . 

While  the  starting  point. for  the  control  system  computations  is 
similar  to  that  for  the  Euler  beam,  our  ^ulalysls  takes  a  somewhat 
different  tack.  We  shall  exploit  the  appearance  of  Hilbert  transfoxms 
in  the  derivation  of  the  spectral  factor  and  the  simple  way  in  which 
these  transforms  can  be  used  to  represent  the  spectral  factor 
appearing  in  the  expression  for  the  feedback  gain.  As  we  shall  see, 
there  are  some  significant  computational  advantages  obtained  in  this 
modification  of  the  Davis-Stenger  procedure. 

In  the  next  subsection  we  describe  the  model  and  compute  the 
transfer  function.  In  the  following  subsection  we  write  the  solution 
for  the  mbsh  dynamical  system  in  terms  of  the  eigenfunctions  of  the 
evolution  operator.  This  provides  an  effective  and  accurate  basis  for 
numerically  simulating  the  (controlled)  system  dynamics.  It  is 
superior  to  numerical  solution  of  the  FDE  model  by  finitq  difference 
methods. 
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4*1  Pynaaies  of  a  Vibrating  Membrane 


The  physical  structure  of  the  prototype  experimental  system  used 
by  Lang  and  Staelin  (I982a,b)  was  described  in  section  1.3*  Ve  shall 
describe  the  mathematical  tmalysis  of  the  system  here.  The  linearised 
equation  describing  the  dynamics  of  the  voltage  controlled  mesh  is 

A  =  -  w  ah  +  ffoZ  h-^  V 

(1)  3t^  “  3a^  **36^  M?  MH^ 

A 

haO  on  boundary  of  [0,L]  x  [O.L] 
h(O,a,0)  »  0,  h^(0,a,6)  «  0 

where  the  boundary  conditions  reflect  the  fact  that  the  mesh  is  pinned 
along  its  boundary.  To  reduce  the  computations,  we  make  a  change  of 
coordinates  to  remove  the  derivative  term  in  (l);  that  is, 

(2)  h  (t)  =*  f  (t)  exp  {•“  t} 

¥e  obtain  3^f  ^  ^  3^f  ^  ^  3^f  ^ 

3t^  M  3a^  M  36^ 

(3)  f=o  on  boundary  of  [0,L]  x  [0,L] 
f  (0,a,B)  -  0  f^  (0,a,B)  =  0 


0  /  W  .  1 

-—2  -“P  ^  ^ 

MH 


q _ t 

.2 


dt"  ®  “  &o* 

(5)  f  a  0  on  boundary  of  [0>L1  x[0,L] 


Qt 


--T  +  y 


*0  at  t 


It  is  possible  to  reduce  the  equation  further  to  the  standard  form  for 
wave  equations;  but  we  shall  not  do  this  since  it  will  complicate  the 
boundary  conditions. 


Taking  the  Fourier  transform  in  (5)t  we  obtain 
Pli).  a.  +  a  ^^-fY^  +  a-u 

(6)  ®  “  Mi  ® 

F=>  0  on  bovindary  of  [0,L]  x  [0,L] 


By  splitting  F  and  U  into  real  and  Imaginary  parts,  respectively,  we 
are  led  to  the  following  equation 


(7) 


a 


3^  H 
3a  ^ 


3^H 


36^ 


+  (  +  w^)  E 


Our  immediate  objective  is  to  find  the  Green's  function  corresponding 
to  the  boundary  value  problem  (7). 


To  accomplish  this,  we  shall 
dimensional  Sturm>Llouville  problem. 


reduce  the  problem  to  a  one 
Consider  the  operator 


(8) 


“  "®a  -  (  Y  +  0)  )  u 

3a^ 

with  u  (0)  =  u  (i)  =  0 


Using  this  in  (7),  we  have 

»8  -I'c.  '■--6“ 

3r 

H  «  0  on  the  boundary 


(9) 


and  80,  the  Green's  function  satisfies  the  FDE 
3^  K 

®8  35*”  -  i-cjK  “  -5  to-  C)  fi  (8  -  n) 

(10) 

/N 

*to,8  boundary  of  [0,i]  x  [0,1] 

Let  ({^(01)  be  the  eigenfunctions  of  the  operator  1^,  Ic  -  1,2,... 

We  shall  compute  the  Green's  function  in  the  form 


(11)  Kto,8,5»n) 


-I 

k«l 


(0)  to) 


We  define  the  weighted  inner  product 
(12)  <  ^  /f  a„  ^  \l)  where  0  =  I0,i]  x  [0,«,1 

Substituting  the  expression  (11 )  for  K  into  (10),  we  obtain 


('3)  "B  ‘■"ic  «>  \  *■'  '  K 


But  ^(a)  are  eigenfunctions  of  L 
(14)  L  <|,j^  (a  )  =  ag^ij^  to) 


Multiplying  both  sides  of  (ijl)  by4i^(a)  ,  k  and  integrating  overoi 


we  obtain 
L 


(15) 


‘  L  “"k  ®’  -  ^0  V  *6  "k  «>  ♦»“> 


to)]  da 


-/q  6  (a  -C  )  6(8  -n  )  do 


which  implies  (using  the  orthogonality  of  4  and  the  properties  of  the 
delta  function) 


(3)  ■  \  ^  ^  ‘  i 


a,t(0)  =  aij^d.)  =0 


This  is  a  classical  Sturm-Liouville  problem. 


The  eigenfunctions  of  the  operator  L  satisfy 


”  V  ° 

It  is  a  simple  calculation  to  shov  that  the  eignevalues  are 
(18)  X  -  ^  1,2,... 

'  t>  aj  «6 

and  the  eigenfunctions  are 


(19)  y  ^  (  n  1^  ot)  n  =  1,2,, 


Now  we  have  to  solve  the  Sturm-Liouville  problem  (16).  Referring 
to  (18),  we  have  to  consider  the  three  cases:  X  >  0,  X  *0,  and 
X  <  0.  Note  that  for  each  fixed  oa  there  is  only  a  finite  number  of 
negative  eigenvalues.  To  solve  (16),  it  suffices  to  solve 

a"  -  Xa  =  -  fi(3-  n) 

(20) 

a  (0)  =  a  («-  )  =0 


We  have  dropped  the  indices  for  convenience  and  the  term  (C)  which 

n 

will  be  handled  by  multiplying  the  resulting  Green's  function  by  the 
same  factor. 


Case  1: 


X-  vr  >  .a  - . 


It  is  a  simple  matter  to  show  that  the  Green's  function 
associated  with  problem  (20)  is 


(21) 


1 

Vish  (U^) 


sh  (  p  (i-n) )  ^sh  (u8)  6<3<ti 

^  sh  4in)  sh  (u  W,"  3) )  T1<6<J, 


where  sh(z)  ■  sinh(z).  To  obtain  the  desired  Green's  function 
associated  with  (16),  we  simply  multiply  (21)  by  <^(0*  gives 

(22)  <3  (8,n)  -  o<6<n 

•  1  p  Sh  (U^l)  )sh  (p  n)  sh(V  (A-3)'n<3<l- 

n  n  C  “  -  n 

where  p  "/X 
n  '  n 

Case  2;  X  -  -P^  <  0 


Arguing  in  a  similar  fashion  we  find  that  the  Green's  function 
associated  with  (16)  for  this  case  is 

sin  fp)  5) 

TTihino” 


shOi  (i-n))  ship  3)  0<3<T1 

n  n 

sh  n)  sh  Oi  (1-3) )  n<6<i 

n  n 


(23) 
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Case  3  s 


A“0 


In  this  ease  the  desired  Green's  function  is 

..  rnr  e 


(24)  G  (6 
n 


[2  sin 


0  (2.-n)  o<$<n 

n  i-B)  n<6<S' 


Note  that  (22)  and  (23)  can  he  given  by  the  same  formula  if  we 
take  y  to  be  a  complex  number  (either  i|v_  I  or  |u  I).  Also,  (24)  can 
be  obtained  as  the  limit  of  either  (22)  or  (23)  0*  However, 
it  is  best  to  split  the  expression  as  above  to  maximize  computational 
efficiency.  The  procedure  we  have  followed  in  representing  the 
Green's  function  has  two  advantages  over  the  classical  expansion  of 
the  Green's  function  in  terms  of  eignenfunctions  of  the  whole  problem 
(as  used  in  Lang  and  Staelin  1982b)  .  First,  it  reduces  the 
expression  of  the  Green's  function  to  a  single  sum  instead  of  a  double 
sum  as  in  the  classical  representation.  Second,  the  series  has  a 
strong  convergence  property.  Since  there  is  only  a  finite  number  of 
negative  eigenvalues,  the  main  part  of  the  series  is  given  by  (22) 
which  can  be  expressed  in  terms  of  exponentials.  This  series 
converges  exponentially  fast.  In  fact,  the  bilinear  eigenfunction 
expansion  for  the  Green's  function  fails  in  this  case  since  *  0  is 
an  eigenvalue  for  an  infinite  number  of  w  's,  and  the  series  diverges 
in  these  cases.  The  complete  expression  for  the  Green's  function  is 


K  fci,B,C,Ti) 
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sin  q-  O 
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nrrC  lurg  ^ 

sin  (  A  )  sin  (  il  )  ji^n)  sh  Oi  3)  CKB<ti 

y  sh(y  t)  . - - *■*> - 

sh  (y  Ti)  sh  (y  (i-B) )  it<B<*- 
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Here  we  have  used  the  notation 


(26) 


Ar!+<*L 


V 


and  tn^l  is  the  integer  strictly  less  thatn^,  including  the  case  when 
is  an  integer.  Also, 


On  not  integer 

(27)  e  In  ]  -  {  ° 

°  ^  integer 

and 


2  2  2  2 
n  TT  Y  +0) 


11*1  f  2  f  • « • 


Finally,  note  that  the  Green’s  function  has  the  symmetry  property 
-  K(C/n,a,6)  for  n<B<Jl 


(29)  Kfci,B,C.n) 


4«2  Transfer  Function  of  the  Controlled  Membrane 

Recall  .the  reduced  form  of  the  model  (7}<  Using  the  definition 
of  the  Green's  function  and  the  symmetry  property  (29)  ve  can  express 
the  transfer  function  of  the  system  as 

(50)  H(a,e,«»))  - //■  fa,6,5,n)  u  (C»n,(D)  dC  dn 

Now  consider  the  piecewise  constant  subdivision  of  the  rectangular 
area  of  the  membrane  as  shown  in  Figure  4«1 


Assume  that  the  contTOl  input  is  constant  over  each  small 


rectangle 

(51)  Vl’  ^ . ** 

That  is, 

(52)  u  (  5,n,(i))  =u^j(u)  for 

Using  this,  the  transfer  function  can  he  expressed 

Hfci,8,c»)  ■  ff  tt  (C,  Tl,  w)  a?  dn 

(33)  N  N 

i*l  j*l  j 

6<rij 

To  compute  the  integral  over  it  is  necessary  to  distinguish 

three  cases !  6  <  n j ,  6  >  ^ ®nd  3e  I  Hj  r  j+^J  • 

f 

Case  1 ;  B  <  n j 


Substituting  the  expression  for  the  Green's  function  (25)  into 

(33),  ve  obtain 
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Bach  of  the  integrals  in  (34)  may  be  computed  in  closed  form.  ¥e 
shall  omit  the  details  and  state  the  final  result 
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Case  2;  6>  n.  , 

3+1 


The  final  result  in  this  case  is 
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Case  3:  ^  6  ^ 
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This  case  may  be  treated  by  substituting  6  for  ri^  in  Case  1  or 
for  Case  2.  The  final  result  is 
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4«3  Solution  of  the  elliptic  system  using  the 
discretized  Green's  function 


We  can  use  the  discretized  Green's  function  as  the  basis  for  an 


effective  algorithm  to  compute  approximate  solutions  to  the  system. 


The  algorithm  is  more  efficient  than  direct  numerical  integration  of 


the  partial  differential  equation.  In  general  terms  the  procedure  is 


as  follows:  Consider  the  complex  elliptic  PDB 


3  H  3  ^  H 

a  — r—  +  b  — —  +  pH  =  -bH  in 


H  =  0  on  aJ2,  n  -  [  0,L]  X  [0  ,L] 


where  a  and  b  are  as  in  (?)  and  p  is  a  complex  parameter  (which  will 
depend  on  the  frequency  O)  ) .  Let 
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^  otherwise 

Assume  that  the  domain  Q  -  [0,l]  x  [o,L]  is  subdivided  into  small 
rectangles  as  in  Figure  4*1 


(40)  =  [x^,  *i+l^ ^j+l^  i=l,...,N, 

Let  h  ••  L/N  and  k  -  L/M.  Assume  that  the  control  is  constant  in  space 
over  the  rectangle  fi..  and  defined  by  its  value  at  the  center.  That 


(41) 


P) 


•’*.  *”7^  . 


i; 


. V*i+i  yj-^^j+i 

u6c,y,p)  =  u^jt>)  =  u  ( - 2 - '  - 2 - ' 


Let  G^j(x,y,p)  be  the  average  Green  function  over 


(42)  G.  .  Oc,y,p)  ■  //o  K  bt»y»C,Tl,p)  d^dn 
.  “ij 


Then  the  solution  to  the  complex  FDB  is  given  by 


N  M 

Hti£»y.p)  -EE  bu  .  te)  G,  .  (>c,y,p) 
i-1  j-1  ^  ^ 


¥e  introduce  the  change  of  notation 


G^j  Cx.y.P)  =  G(x,x^,x^^^,y,y^,yj_j^j,p) 


which  will  be  useful  in  programming  the  algorithm.  The  function 
Gij(x,y,p)  will  be  given  by  the  following  expressions: 


A  ^ 


(45a) 


G  (x.yrp)  =  -  (Gj^  Wjf  p)+G2^"J^»  y  <  yj 


(45b)  =  G^(x,x.,x.^^,y,  y^^yj+i^p)  ^G^  (-)  ,  y .  <  y 


V  '1'^^  V?^?i^?.vV-V.V?.*.V//.  "Vs*.  *\  -VAV.  A  A  A  •*.  •*.  •••  •“.  •  .  A  A  A  A  ^ 


(45c) 


G 


L-y,p) 


J.XJ. 


C>c,y,p)  =  G  C>c,x^x^^j^,y,y^^y,p)  +G  tK,x^,x^^j^,L-y,  Wj+i. 


Where 


Q.  (-)  =  {  ?  1  - 2^ - sr- 

nTTV^shOxL). 


sin  {^)  sh  tU^(L-y)]  } 


(46a) 


nwx.  ,  1 

X  {cos  ^ -  cos  ( — — — )}  x  ^ch  (lij^yj^j^)-ch  (y^y^)} 


(46b) 


2  “o"'* 

^  ^  U 

x-|  (  -  y^  ) 

2  ^j+i 


CL-y)  •  [cos  t 


nQTT 


n  IT 

X.)  -cos  {“7 —  *-• 

1  li 


i+1 


)] 


The  conputatiooal  al^oritho  based  oa  this  representation  of  the 
solution  is  given  in  Figure  4<2. 


4>4  Eigenfunction  expansion  of  the  solution 
of  the  generalized  wave  equation 


In  this  subsection  we  shall  use  an  eigenfunction  expansion  to 
solve  the  system  (1)  as  ,  an  evolution  system  in  time.  Using  the 
exponential  transformation  (2)  we  can  rewrite  the  system  as 


-^-2  -  a  +b  +dg  +bu  in  fZ 

3t  ax'^  3y^ 

(47) 

g»0  on  an,  g(x,y,tQ)-  g^  (x,y) ,  g^(x,y,t^j)=  Cx,y) 
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V  vv^.* 


i 


Here  we  use  general  boundary  conditions  since  we  want  to  use  the 


solution  to  simulate  the  response  of  the  membrane  system  to  inputs 

based  on  its  state  at  a  given  time.  Our  solution  is  expressed  in 

teins  of  an  eigenfunction  expansion.  Let  (x,y)  be  the 

mn 

eigenfunctions  of  the  operator  in  (47)  and  suppose  that  the  data  and 
the  control  function  have  the  form 


u  Oc,y,t)  »  ^  u  (t)  0  bt,y) 
ia#n  mn  mn 


go  -  2  (K.y)  g  &c»y,tQ)  -  z  t>c,y) 

niyn  X  mpn 


Ve  look  for  the  solution  in  the  form 


(49)  m,n 


This  leads  to  an  eigenvalue  problem 


a  +  <|(K  mX  i 

,  .  -  2  ^2  mn  ^iim 

(50)  3x  3y 

(p  <>-  0  on 
mn 


and  an  initial  value  problem 


a  (t)  »  X  a  (t)  +b  u  (t) 
mn  mn  mn  mn 

“mn  "  >^mn  '  “mn  “^mn 

mn  mn  mn  mn 


A  simple  calculation  determines  the  eigenfunctions 


(52) 


corresponding  to  the  eigenvalues 


(53)  -  1.  , 


The  initial  value  problem  has  different  solution  forms  depending 


on  the  sign  of  X  .  Let  t_  be  the  initial  time  and  let  t.  -  t.  T, 
mn  0  £0 


where  T  is  the  sampling  period.  VhenX  <  0,  we  have 

BBl 


nm  h 

a  (t  _)  ■  r  cos  <li  T)  +  -  sin  IU  T)  +  - —  •  u 

nn  f  mn  mn  y  mn  y  mn 

mn  mn 


l-cos  y 


rm 


(54a)  .  sin  u  T 

a  (t_)  «  -  r  y  sin  (y  T)  +  V  cos  T)  +bu  mn 

mn  f  mn  mn  mn  mn  mn  mu 


mn  mn 


mn 


mn 


mn 


'^mn 


■  0,  we  have 
mn  ' 


o  (t-)-  r  +  v__T  +bu  ~ 
nn  f  mn  mn  mn  2 


(54b) 


ct  (t^)  -  V  +  bu  T 
mn  £  mn  mn 


When  X  >  0,  we  have 
mn 


a__  (t*)  =  ch  sh  «nn  eh 


-mn  "-f ' 


mn  mn 


mn 


mn 


nn 


nn 


(54c) 


a  (t.)  *y  r  sh  (y  T)  V  ch  (U  T)  +bu - 

mn  f  mn  mn  mn  mn  mn  mn  u 

non 


These  expressions  permit  us  to  solve  the  membrane  system  from  any 
initial  conditions  for  any  control  input  satisfying  the  condition  that 
it  be  constant  over  the  area  of  the  mesh  element  in  Figure  4.1. 


4o  Spectral  factorization  using  the  Hilbert  transform 


Computation  of  the  optimal  stabilizing  gain  based  on  the  Davis  > 
Stenger  spectral  factorization  procedure  used  for  the  Euler  beam 
problem  proved  to  be  infeasible  in  treating  the  tvo  dimensional 
system.  The  computational  effort  was  too  high  to  be  practical;  and 
it  was  necessary  tc  redesign  the  algorithm  to  achieve  greater 
numerical  efficiency.  The  key  idea  was  to  use  the  Hilbert  transform 
to  represent  the  spectral  factors  arising  in  the  calculation.  This 
permitted  the  use  of  fast  and  efficient  numerical  Fourier  transform 
techniques  in  the  gain  computation  algorithm. 


Given  a  Fourier  pair  (f,F)  in  xL^ 


C55)  f  (x)  -  ^  C  P(8)ds,  P(s)  -  -I-  ^e^**f(x)dx 


the  Hilbert  transform  is  defined  by 


m  (t) 


lfc)_  a,  51,2 
TT  -»  t-r 


where  the  Integral  is  taken  in  the  Cauchy  principal  value  sense. 
Giv.en  a  complex  function 


(57)  T(s)  -  R(s)  +  iF(a) 


(58)  HT(  oj)  -  HB(  u)  +  iHF(  w) 


f 

>  j" 

>  > 

■r 


“V  -2*- I 


The  inverse  Fourier  trsnafora  of  HT  satisfies 
(59)  F"^  (ht)  -  -isgn(t)  F"^  {t} 


If  we  define 


4  CP-iHP),  T  »  4  (  14-iHC) 

(60)  +2  -  2 

then 

t>0 
t<C 


r 


(61) 


t>o 
C)  t<  0 


,  p‘^«r  )  - 


^P"^  CT) 


This  means  that  T^  (respectively  T  )  is  the  projection  of  T  onto  the 
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space  (respectively  L  )  with  the  obvious  intepretation  of  and 
2  * 

L  . 


Now  considor  a  function  f(x)  satisfying  the  condition 


(62)  lim  f  (x)  «  1 

I  X  !  "  2 

if  ^  (z)  >•  f(z)  -  1,  then  e  L  .  Applying  the  previous  results  to 
the  function 


(63)  -  log  [f  be)] 

(with  the  proper  choice  for  the  branch  of  log),  then 

(64)  f(x)  -  f+(x)*f-(x) 


where  f*  and  f~  are  the  causal  and  anti-causal  spectral  factors. 


respcctivaly.  For  Instanc*, 

(65)  f“(x)  -  f(x)  expi-i  Hlog  f(x)} 

2 

This  sxprsssion  Is  ussd  directly  in  the  computation  of  the  optimal 
stabilising  feedback  gain  by  spectral  factorization. 

4*6  Gain  computations 

In  abstract  form  the  control  system  (1 )  takes  the  form 

X  (t)  ■  Ax  +  Bu 

(66)  y  (t)  -  Cx  ,  X  (0)-  Xq 

The  transfer  function  is 
(5Y)  GCiw)  »  CRttujA)  B 

where  R(iu)  ;A)  is  the  resolvent  operator 

(68)  R(iw;A)  -  (iWI-A)"^ 

If  we  compute  the  spectral  factor  of  I  *  G*G 

(69)  f'  (ii*))  F*  (iW)-  tt+G*G)  Uu) 

then  the  gain  is  given  by 

[B*K]  (Xq)  -  G*ttW)CRaw,A)  X 


(70) 


In  the  present  case  the  components  are 


And  ve  shall  take 


(72)  o 


^2j 


^  tt.^  a  L^) 


which  allows  for  weights  on  both  the  displacement  and  the  velocity 


The  resolvent  operator  is  computed  by  solving  the  system 


(73) 

(:)■ 

which 

leads 

to  the  system 

3^f 

3^f 

+  b  — =■  -  (s 

3y 

(74) 

f=0  on  3fl,  g"sf-u 


The  associated  transfer  function  is 


aV  2 

a  •«>  — ~  +  «i)  +ic  )  g,  —1 

(76)  9*2  3y2 

gj^  on  3fi  ,  g2  -i  g^ 

Therefore,  if  ve  define  H(ia))  ■  (I  Cr6)(ia»  )  we  obtain 

(77)  H(ia))  -  1  ♦  (K^  +u»^K^)  g^(j  u) 


Let  the  state  of  the  system  be  defined  by 


(78) 

Let 

[:]-  [:] 


and  let  satisfy  the  eonplez  PQE  (76).  Then  the  stabilizing 
feedback  gains  are  given  by 


B*kL  I-  ^  f"  ttu)  ‘^g*  (iw)K(iu)  (kJ  + 


+  ia>K^  F  Ciw) 


4.7  Software  development  and  control  system  performance 


Fortran  code  has  been  developed  to  implement  the  control  gain 
computation  (80)  using  the  Hilbert  transform  representation  (65)  for 
the  anti-causal  spectral  factor  of  the  return  matrix.  The  code 
includes  the  simulation  algorithm  for  the  system  response  to  the 
control  as  shown  in  Figure  4.2.  Testing  of  the  code  has  been  carried 


•*! 


out  for  a  few  values  of  the  weights  in  (72).  The  results  indicate 
that  the  stabilisation  of  the  systea  is  significantly  enhanced  by  the 
control  action.  Since  we  take  the  damping  in  (1)  to  be  positive,  the 
system  is  stable  for  small  enough  bias  voltages.  Modest  values  of  the 
control  weights  (10  on  a  normalised  scale),  produce  a  20%  improvement 

f 

in  the  settling  time  of  the  system.  Extensive  testing  will  require  an 
upgrading  of  the  numerical  routines  in  the  code  to  produce  faster 
solutions  of  the  system. 


Figure  4.2  Algorithm  for  control  gain  confutation  emd  sinrulation 


5.  Homogenization  of  Regular  Structures 

It  is  now  generally  accepted  that  large,  low  mass  lattice 
structures,  e.g.,  trusses,  are  natural  for  space  applications.  Their 
large  si^e  and'  repetitive  infrastructure  require  special  techniques 
for  structural  analysis  to  cope  with  the  large  number  of  degrees  of 
freedom.  As  Noor,  Anderson,  and  Greene  (1978)  point  out,  continuum 
models  provide  a  simple  means  for  comparing  structural  characteristics 
of  lattices  with  different  configurations,  and  they  are  effective  in 
representing  macroscopic  vibrational  modes  and  structural  response  due 
to  temperature  and  load  inputs.  Our  approach  to  the  construction  of 
such  models  is  presented  in  this  section.  In  the  next  two  sections  we 
consider  the  problems  of  control  and  state  estimation  in  combination 
with  the  construction  of  continuum  models.  Ve  shall  begin  with  a  few 
remarks  on  related  work  on  continuum  models  in  the  recent  structural 
mechcmlcs  literature. 

Noor,  et.  al.  (1978)  use  an  energy  method  to  derive  a  continuum 
approximation  for  trusses  with  triangular  cross  sections  in  which  the 
modal  displacements  of  the  truss  are  related  to  a  linearly  varying 
displacement  field  for  an  equivalent  bar.  In  (Dean  and  Tauber  1959) 
and  (Renton  1969),  exact  analytical  expressions  for  the  solutions  of 
trusses  under  load  were  derived  using  finite  difference  calculus.  ^ 
expressing  the  difference  operators  in  terms  of  Taylor's  series  Renton 
(1970)  was  able  to  derive  continuum  approximations  to  the  finite 
dlfferer>ce  equations  resulting  in  expressions  for  equivalent  plate 
stiffnesses,  for  example.  In  a  recent  paper  Renton  (1984)  used  this 


approach  to  give  equivalent  beam  properties  for  trusses,  and  this 
complements  the  earlier  work  of  Noor,  Anderson  and  Greene  (1978)  and 
Nayfeh  and  Hefzy  (1978).  (See  also  (Anderson  1981).) 

In  these  papers  a  continuum  model  is  associated  with  the  original 
(lattice)  structure  by  averaging  the  parameters  of  the  lattice  over 
some  natural  volume  (e.g.  of  a  "cell"  of  the  structure)  and 
identifying  the  averaged  parameter  value  (mass  density,  stress  tensor, 
etc.)  with  the  corresponding  distributed  parameter  in  the  continuum 
model.  A  specific  form  for  the  continuum  model  is  postulated  at  the 
outset  of  the  analysis;  e.g.,  a  truss  with  lattice  structure  will  be 
approximated  by  a  beam,  with  the  beam  dynamical  representation  assumed 
in  advance.  While  this  approach  has  an  appealing  directness  and 
simplicity,  it  has  some  problems. 

First,  it  is  very  easy  to  construct  an  example  in  which  the 
"approximate  model"  obtained  by  averaging  the  parameters  over  a  cell 
is  not  a  correct  approximation  to  the  system  behavior.  This  is  done 
in  subsection  9.1.  Second,  the  averaging  method  (averaging  the 
parameters  over  space)  does  not  apply  in  a  straightforward  way  to 
systems  with  a  random  structure,  since  the  appropriate  averaging 
procedure  in  this  case  may  not  be  obvious.  Third,  the  method  cannot 
be  naturally  Imbedded  in  an  optimization  procedure;  and  controls  and 
state  estimates  based  on  the  averaged  model  may  not  be  accurate 
reflections  of  controls  and  state  estimates  derived  in  the  course  of  a 
unified  optimization  -  averaging  procedure.  The  method  does  not 
provide  a  systematic  way  of  estimating  the  degree  of  suboptimality  of 
controls  and  state  estimates  computed  from  the  idealized  model. 


In  this  work  we  use  a  totally  different  technique  called 
homogenization  from  the  mathematical  theory  of  asymptotic  analysis  to 
approximate  the  dyneunics  of  structures  with  a  repeating  cellular 


structure.  Homogenization  produces  the  distributed  model  as  a 
consequence  of  an  asymptotic  analysis  carried  out  on  a  rescaled 
version  of  the  physical  system  model. 

Unlike  the  averaging  method,  homogenization  can  be  used  in 
combination  with  optimization  procedures;  and  it  can  yield  systematic 
estimates  for  the  degree  of  suboptimality  of  controls  and  estimators 
derived  from  idealized  models.  Results  to  this  effect  are  given  in 
sections  6  and  ?•  While  our  results  are  preliminary,  they 
nevertheless  demonstrate  the  feasibility  of  the  method;  and  they 
suggest  its  potential  in  the  analysis  of  structures  of  realistic 
complexity. 

In  this  section,  we  first  give  an  example  illustrating  some  of 
the  subtleties  of  homogenization;  then  we  discuss  homogenization  for 
abstract  hyperbolic  systems;  then  we  illustrate  the  applications  of 
homogenization  theory  by  deriving  a  diffusion  approximation  for  the 
thermal  conductivity  of  a  (random)  lattice  structure.  In  the  final 
example  we  derive  a  homogenized  representation  for  the  dynamics  of  a 
lattice  structure  undergoing  transverse  deflections.  We  show  that  the 
behavior  of  the  lattice  is  well  approximated  by  the  Timenshenko  beam 
equation;  and  we  show  that  this  equation  arises  naturally  as  the 
limit  of  the  lattice  dynamics  when  the  density  of  the  lattice 
structure  goes  to  infinity  in  a  well  defined  way.  The  mathematical 
analysis  used  in  the  derivations  is  based  on  the  book  (Bensoussan, 
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Lions,  and  Papanicolaou  1978)  and  the  paper  (Kunnemann  1963).* 
3»1  A  one-dioensional  example 


From  (Bensoussan,  Lions  and  Papanicolaou  1978)  we  have  the 
following  example: 

,  -k 

(I) 

(jCq)  -  0  =  u^  hc^) 


where  a(y)  is  periodic  with  period  y  ,  a(y)  >  a  >  0,  and  a^(x)  ■ 

0  ~ 

Siij/e)  •  It  is  simple  to  show  that 


(2)  llu^  11^  A  I  S 

and  so,  u^  -»-u  weakly  in  the  space  H*" .  Moreover, 

(3)  a^'-^MCa)  ‘iy 

^0  0 

and  it  is  natural  to  suppose  that  u^->  u  with  the  limit  defined  by 


_  _f _ [  M(a)  --u  (x)  ]  -  f(x)  X  e  CXq,Xj^) 

“  dx  “ 

u  (X  Q )  0  ■  u  (x^  )  ’ 

This  is  untrue  in  general  (Bensoussan,  Lions  and  Papanicolaou 


*¥e  are  grateful  to  Professor  George  Papanicolaou  for  bringing 


Kunnemann 's  paper  to  our  attention 


1978,  pp«  8-10) •  Th«  correct  limit  is  given  by 


(5)  '  5?  '*  ‘*1> 


with 


(6) 


In  general,  M(a}  >  a;  and  so,  the  error  is  identifying  the  limit,  (4) 
versus  (5),  is  fundamental. 


The  system  (4)  corresponds  to  averaging  the  parameter  a  (z)  over 
a  natural  cell;  a  procedure  similar  to  that  used  in  (Noor,  Anderson 
and  Greene  1978),  (Nayfeh  and  Hefzy  1978)  and  (Aswanl  1982)  to  define 
continuum  models  for  lattice  structures.  As  (9)  shows,  the  actual 
averaging  process  can  be  more  subtle  than  one  might  expect,  even  for 
simple  problems. 


To  see  how  (9)  arises,  let 


(7)  5^(x)-a'^(x) 


£  2 

Then  5  (x)  is  bounded  in  L  (Xq,Xj^)  and  it  satisfies 


(8)  -  h  W-f  W,  xe 

091 


One  can  show  that  V"  (i)  has  a  strong  limit C  (x)  in  J?  j  so 


(9)  « 


weakly  in  L^(x  ,x  ).  But 


.  jr.  >-fm  ■-.-  •».-  • 


*•  . 
t 


pt ,  e  du 

(10)  «  /•  -  al 


SO 


(11) 


and  since 


(12)  -i-£ 


we  have 


^  r/Mr—  1 

(13)  “  dx  ^  a  dx  ^ 


N.'. 


and  the  limit  u  is  weak  in  H^(i^tX2)  (the  Hilhert  space  with  norm 

defined  hy  (2)). 


c: 


i 


This  example  illustrates  the  pitfalls  associated  with  simplistic 
averaging  procedures. 

5*2  Homogenization  of  Wave  Equations 

As  shown  in  Part  I  of  this  report,  hyperbolic  (wave)  equations 
are  the  most  natural  models  for  the  dynamics  of  flexible  structures. 
General  techniques  for  homogenization  of  wave  equations  are  available 
(Bensoussan,  Lions  and  Papanicolaou  1973).  The  precise  form  of  the 
homogenization  procedure  depends  on  the  scaling  of  the  physical  model; 
l.e.,  the  dependence  of  the  system  characteristic  features  on  small 
parameters.  Since  this  is  a  sensitive  modeling  issue,  we  shall 
discuss  it  in  some  detail.  It  may  be  necessary  to  consider  several 
scalings  to  determine  the  one  most  suitable  for  a  given  class  of 


•:nVA>wLv’ 
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flexible  structures* 


Consider  the  Kleln-Gordon  equation 


u^^{t,x)-V[c  (x)  7u(t,x)l  -Wbc)  u(t,x) 


u(0,x)«f  (x),  (0,t)-g  (x)  ,  t>0,  xe  R 

where  f  and  g  are  smooth  functions  of  compact  support,  c  (z)  >  0  and 
V(z)  ^0  are  smooth.  Suppose  c(z),  W(z),  f(z),  and  g(z)  depend  on  a 
small  parameter  e  >  0.  For  example,  suppose  c  and  ¥  vary  slowly  with 
X,  so  we  have 


(15)  c^«c^  (£  X),  W-W(e  X) 


and  suppose 


(16)  f  -  f®,  g  -  g' 


Let  u  be  the  solution  of  (14)  with  (15)  (16).  The  behavior  of  u  as 
^ ’t’ 0  In  (14)-(16)  Is  trivial  If  we  consider  (t,x)  fixed.  However,  If 
(t,x)  become  large  as  e  -^0,  then  an  Interesting  limiting  behavior 
emerges  (Bensoussan,  Lions  and  Papanicolaou  1978,  Chapter  4). 


To  see  this  most  clearly.  It  Is  necessary  to  rescale  t  and  x  as 


(17)  t'  =  et , 


Dropping  the  primes,  this  yields 


(t,x)=7*[c^  (x)7u^  (t,x)]-  ^  w  (x)  (t,x) 

(18)  ^ 

U  (0,x)«f  ,  u^  (  ,x)^  X  £  R*',  t  >0 
Notice  that  In  rescaling  a  system  with  slowly  varying  coefficients. 


one  produces  a  systma  trlth  a  large  coefficients 


To  complete  the  specification  of  the  problem  it  is  necessary  to 
define  the  dependence  of  the  data  f  ^  and  g  ^  on  e .  The  various  choices 
strongly  Influence  the  final  form  of  the  limiting  system.  Follovlng 
(Bensoussan  et  al.  1978),  we  shall  distinguish  three  classes  of  data. 

Case  1 :  Low  frequency  problem 

In  this  case  and  g^  have  asymptotic  power  series  expansions 

(x)  f  be)  +  E  f.  be)  +••• 

(19)  ° 

bt)  gQ  bt)  +  c  be)  +... 

where  the  terms  in  the  expansion  are  smooth  functions.  To  produce  a 

problem  with  finite  energy  as  e  ***0,  it  is  necessary  to  take  TqCx)  -  0 

in  (19).  (See  (Bensoussan  et  al.  1978  p.  941)*)  Since  the  problem 

2 

is  linear,  it  is  not  necessary  to  have  the  expansions  begin  with  t  . 

The  analysis  of  this  problem  is  comparatively  simple,  and  the 
limiting  behavior  as  e  -*-0  is  elementary. 


Case  2:  High  frequency  wave  propagation  in  a  slowly  varying  medium 

In  this  case  the  data  take  the  apparently  specialised  form 

bt)  exp  0  i  S  bt)  /€]  <x) 

bej  -  exp  [  is(x)/el  g^  (x) 
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where  S(x)  is  real- valued  and  saooth,  and  and  gC  are  complex-valued 
and  smdoth  and  have  asymptotic  expansions  like  (t9)>  Vote  that  the 
data  in  (20)  are  complex- valued.  Since  the  problem  (18)  is  linear, 
both  the  real  and  imaginary  parts  of  u  are  solutions. 


To  understand  the  physical  significance  of  the  second  case, 
suppose  the  "phase  function"  S(x)  *■  k  x,  where  k  is  a  constant  vector 
and  dot  stands  for  the  inner  product  in  R^.  In  this  instance  the  data 
in  (20)  are  spatially  modulated  plane  waves  with  rapidly  varying 
phase.  As  shown  in  (Bensoussan  et  al.  1978)  virtually  all  cases  of 
interest  (different  scalings  leading  to  nontrivial  limiting  behavior) 
can  be  analyzed  in  terms  of  this  case.  For  instance  it  is  possible  to 
treat  the  case  of  wave  propagation  in  slowly  varying  media  with 
spatially  localized  data  or  the  form 


(21) 


l-n/2 

Oc)  -  e 


e 

g 


(x)  «  e 


-n/2 


f(>c,|) 

f  ) 


where  f  and  g  are  smooth  functions  of  compact  support  in  x  cmd  y 
(  x/e),  and  (x,y)  eR^  .  The  scaling  on  the  right  in  (21)  is  chosen 
so  that  the  terms  are  of  order  one  as  0.  Problems  with  forcing 
functions  and/or  inhomogenous  boundary  conditions  can  also  be  treated 
by  essentially  the  same  method. 


The  treatment  of  case  2  given  -  (Bensoussan  et  al.  1978)  is 
based  on  the  ideas  of  geometric  opt^  ,3.  The  solution  u^  is  sought  in 


(22) 


u  (jc»t)  -  e  e 


is(x,t)/e^e 


V  (x,t)  -Vq  (x,t)  +  ev^  (x,t)  +  ... 

Inserting  (22)  into  (18)  and  equating  coefficients  of  equal  powers  of 
e  leads  to 


(23)  »1  ’'l  *  *2  ''O  -  “ 


Aj  Vj  +  Aj  .  Aj  »o  -  0 


where  o  22 

-  -c^^  (  7  S  +  (S^)^  -  N 


(24) 


A2  -  -  2  i  St  3^  V  •  (c^  7  s  )  +  ic^7  S  •  7  -  i  Stt 


A-  -  7  •  fc*'  7)  -  3, 


The  analysis  of  these  equations  leads  to  the  reduced  model  for  the 
system  behavior. 

From  the  first  expression  in  (23)  we  see  that  for  v^  not  to  be 
identically  zero  we  must  have 

(25)  St-  [  CV  S  )^  +  W]  -  0 

This  is  the  eikonal  (or  Hamilton-Jacobi)  equation,  a  nonlinear  first 
order  PDB  which  controls  the  evolution  of  the  phase  function.  It  may 
be  solved  in  terms  of  a  system  of  nonlinear  ODE's  (Hamilton's 
equations)  for  the  "rays"  cuid  ’’momenta"  associated  with  the 
propagation  of  energy  by  the  system.  (See  (Bensoussan  et  al.  1978, 
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The  choice  (25)  for  the  phase  makes  the  operator  Identically 
zero.  Using  this  in  the  second  equation  in  (  3  )  leads  to  the 
"transport  equation" 

(26)  ^  ^0  "  0 

Its  analysis  leads  to  expressions  for  the  propagation  of  energy  in  the 
system. 

The  case  of  spatially  localized  data  (21 )  may  be  treated  using 
the  same  techniques  in  combination  with  the  method  of  multiple  scales. 
Ve  shall  not  develop  the  analysis  of  these  general  systems  in  more 
detail.  Rather,  we  shall  turn  our  attention  to  the  construction  of 
continuum  models  for  lattice  structures. 

5.5  Continuum  approximations  for  lattice  structures 

In  this  section  we  shall  apply  homogenization  and  the  associated 
asymptotic  emalysls  to  derive  continuum  approximations  for  two 
different  types  of  problems.  In  the  first  case  we  show  that  the 
problem  of  thermal  energy  conduction  in  a  lattice  can  be  well 
approximated  by  a  diffusion  process  in  the  macroscopic  scale.  In  the 
second  case  we  show  that  the  in- plane  macroscopic  (2  dimensional) 
motions  of  a  simplified  truss  model  can  be  well  approximated  by  the 
Timoshenko  beam  system. 


3>3*1  Effective  conductivity  of  a  periodic  lattice 


A.  Problem  definition 

Let  Z  ■  {0,^1  ,^2, . . . }  and  2?  *  Z  x.  ..x  Z  (d  times)  be  a 

d-dimensional  lattice.  Let  e  >0  be  a  number  small  relative  to  1 .  We 
want  to  describe  the  effective  conduction  of  thermal  energy  on  the 
e -spaced  lattices  Let  e^  -  (0,...0, 1,0, ...O)^  ifith  1  in  the  ith 

position,  1  -  1,2,...,d.  If  x  is  a  point  in  Z*^,  then  x  ^  £ej^, 

1  ^  if.  d,  are  the  nearest  neighbors  of  x.  Let  a^^(x),  x£  Z*^, 

1  .<  i  .<  d,  be  the  two  functions  defined  on  the  lattice,  and  assume 

(26a)  (x)  Cx)=a^- (x+e^) ,  x  e  Z^,  1  <  i  <  d 

(26b)  0<A  i  a^(x)-  B<»xe  Z*^,  1  »  i  -  d 

(26c)  ,a^  (x)  is  periodic  with  period  ^  1  in  each  direction,* 

1  <  i  <  d. 


Next  let 

(27)  (x)  ■  (x/e) ,  xCeZ**,  1  S  i  S  d 

Equation  (26b)  means  that  the  conduction  process  is  reversible 
and  that  the  conductivity  a^  (x)  is  a  "bond  conductivity",  i.e., 
independent  of  the  direction  in  which  the  bond  (x,x  e^)  is  used  by 
the  process.’  Equation  (27)  means  that  the  configuration  of  bond 
conductivities  one, ^  is  simply  njL^.(’)  on  Z*^  "viewed  from  a 

distance."  Assumption  (26c)  imposes  a  regularity  condition  on  the 
le  period  may  be  different  in  different  directions. 


K,*. 


pbysics  of  the  conduction  process.  An  assumption  like  this  is 
essential  for  existence  of  a  limit  as  e  -^0.  In  one  dimension  the 
situation  is  illustrated  in  Figure  5.1.  A  system  similar  to  this  was 
treated  by  Kunnemann  (1983)  with  random  bond  conductivities. 
Brgodicity  replaced  periodicity  in  (Kunnemann  1983). 


-2t  -  « 


t  11  3<  ff  . 


Figure  S.l.b.  Conductivity  on  c-scaled  lattice,  y  -  ex,  with 

X  e  Z  and  period  e£. «  6  . 


One  can  associate  with  this  system  a  random  (jump)  process 

{X^(t,x),  t  0,  x£  CZ^}  on  the  e-spaced  lattice*.  In  effect,  as 

c  0,  {X^}  converges  to  a  Brownian  motion  on  the  lattice;  and  the 

main  result  of  the  analysis  is  an  expression  for  the  diffusion  matrix 

*Deflnltlon  of  this  process  is  not  necessary  for  the  analysis,  but  it 
bolsters  the  intuition. 


>!.'■  .N'. 


Q:  ■  of  this  process.  This  matrix  describes 

the  macroscopic  diffusion  of  thermal  energy  in  the  system.  It  is  the 
effective  conductivity. 

¥e  shall  carry  out  the  asymptotic  analysis  of  this  system  in  the 
limit  as  e->0  using  the  theory  of  homogenization.  Let 

(V^®'“u)  (x) ^  [u(x-ce^)-u(x)  1 

(28)  (x) i  lu (x+£e^) -u (x) 1 

x€eZ^,  1  i  i  i  d, 

for  any  u  square  summable  on  ^  or  square  integrable  on  with  e^ 

the  ith  natural  basis  vector  in  R^.  Then 

3u^(t,x)_  _  ^  c- .  ^  e+^e j 

(29)  i-1  ^  " 

s-L'^u®(t,x) 

is  the  diffusion  equation  on  the  e-spaced  lattice  with  density  u(x) 
and  conductivity  a^.(z/£)«  Our  construction  of  an  effective  parameter 
representation  of  the  thermal'  conduction  process  as  e  0  will  be 
based  on  an  asymptotic  €malysis  of  (29)  using  the  methods  in 
(Bensoussan  et  al.  1978). 

Remark:  Although  we  shall  not  use  probablistic  methods  in  the 

analysis,  the  associated  probabilistic  problem  has  a  great  deal  of 

£ 

intuitive  appeal.  The  operator  L  may  be  identified  as  the 
infinitesimal  generator  of  a  pure  jump  process  X^(s)  in  the  "slow" 
■  e  t;  cf.  (Breiman  1968).  Moreover,  L  is 


time  scale  s: 


'  .  •  1  •  r  .  >1 .V-  > 


d 

selfadjoiut  on  Z  with  the  inner  product 
(50)  “'«>• 

Hence,  the  backwards  and  forwards  equations  for  the  process  X  (s;  are, 
respectively, 

-  tLS^(y.t[*)l  (X) 

O't 

(31) 

-  [LS^(-»t|x)]  (y) 

So  the  process  is  "symmetric"  in  the  sense  of  Harkov  processes 
(Breiman  1968). 

The  asymptotic  analysis  of  (29) »  when  interpreted  in  this 
context,  means  that  as  the  bond  lattice  is  contracted  by  e  and  time  is 
sped  up  by  e  ,  the  Jump  process  (X  (s)}  approaches  a  diffusion 
process  with  diffusion  matrix  Q.  In  other  words,  on  the  microscopic 
scale  thermal  energy  is  transmitted  through  the  lattice  by  a  Jump 
process;  but  when  viewed  on  a  macroscopic  scale  the  energy  appears  to 
diffuse  throughout  the  lattice.  The  microscopic  physics  are  described 
in  (Kirkpatrick  1973)  and  (Kittel  1976). 

Because  the  basic  problem  (29)  is  "parabolic",  we  can  introduce 
the  probabilistic  mecheuiism  and  make  use  of  it  in  the  analysis.  In 
the  "hyperbolic"  problems  we  treat  later,  this  device  is  not 


available 


B.  Asyn^totic  analysls-homgenization 


The  essential  mathematical  step  is  to  show  strong  convergence  of 
e 

the  semigroup  of  L  ,  say 


(32) 


(t) : *exp (L^t) 

-»•  T(t)  ;=exp(Lt) 

e+0 


and  identify  the  limiting  operator 


(33)  L=  E  q.  .  ^ 


i.j=l 


This  is  accomplished  by  proving  convergence  of  the  resolvents 


(34)  for  a  >  0,  [ -l^  +a]  [-L+a]’^ 


That  is,  if  f  is  a  given  function  and 
u^(-) :=[-L^+a]~^f 

(35) 

U('  )  :  =  f-L+a]“^f 
then  u  u  ( in  some  sense ) . 

The  method  of  multiple  scales  developed  in  (Bensoussan  et  al. 

1978)  will  be  used  to  prove  the  limit.  Because  the  conductivities 

a  (x)  in  (29)  do  not  depend  on  time,  we  may  work  directly  with  L*" 
i 

rather  than  the  parabolic  PDE  (3)  (of.  (Bensoussan  et  al.  1978 
Remark  1.6,  p.  242.).  The  method  of  multiple  scales  is  convenient 
because  it  is  a  systematic  way  of  arriving  at  the  "right  answers"  - 
something  which  is  not  always  simple  in  this  analysis. 


Following  (Benaousaan  et  al.  1978)  and  bearing  in  mind  (35) t  we 


conaider 

(36)  (l.V)(x)-f(x) 

e 

with  boundary  conditiona  to  be  apecified  later.  Ve  ahall  look  for  u 
in  the  form 

(37)  u^(x)-UQ(x,|)+eu^  (x,|)+€^u2(x,|)+. . . 

d 

with  the  functiona  u^(x,y)  periodic  in  y  ee  Z  for  every  J  -  0,1,.... 
(Aa  it  tuma  out  the  boundary  conditiona  are  aomewhat  irrelevant  to 
the  conatruction  of  "l^lght  anavera.")  To  preaent  the  computationa  in  a 
simple  form,  it  la  convenient  to  Introduce  y  -  x/e,  to  treat  x  and  y 
as  independent  variables,  and'  to  replace  y  by  x/e  at  the  end. 

g+ 

Recall  the  operators  V  from  (29)*  Applied  to  a  smooth 
function  u  >■  u(x,x/e)f  we  have 


(x,y)  *  i  [u(x-ee^,  y”ej^)-u(x,y)] 

»  —  [u  (x,y-e^)-u  (x,y)) 

(38) 

+  i(u(x-  e.,  y-e.)  -uCx.y-e.)] 
e 


2 

i  (  V"  u)  Cx,y)-  1^  Cx,y-e^)  +e  |  |^2  ’ 


where  on  functions  ♦(y) 


(39)  (^)  fy)  m  ^  (y) 
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Defining 


(40)  )(y)  -  (j><y+eiJ  - 


we  also  have 


(41) 


(x.y) 


1  (7*!  u)  (x,y)+  (x,y+«^)  +ei 

£  1  " 


Now  we  substitute  (37)  into  (36)  and  use  the  rules  (39)  (40). 

Equating  coefficients  of  like  powers  of  e,  this  leads  to  a  sequence  of 

equations  for  u  ,  u  Specifically,  (using  the  summation 

0  1 

convention) 


(42) 


(L^u^)  (x,y)  -  (y)  ] 

*  ’  “2  t  ^i^y)  7^  Ujjbc.y)] 

“7^  I  (y)  —  c*»y+«£>  1 

•J  e  7  ®  [a^  (y)  — -2-  6c,y+e  )  +0  (e) 
3x^  ^ 

-  E  ^i  I  (y)  7^  u^  tx,y)  J 
£- 

air  +o(e;- 

1 

-7i  [a^(y)  7^  u^  (x,y)J  +0  (e)  -  f  (x) 


That  is,  labelling  each  term  by  its  order  in  e 

(43)  ('='^)  »o'  ■  ° 

/  \  •  3u^  _  M 

(44)  {e"^)  e7^^‘  [Si  (y)  ^  (x,y+e^)]+7.  Ia.(y)7.  u^(x,y)]-  0 

and  (recall  e7^'  -  0(1 )  in  e  ) 


2^^  6c,y+e^)] 

0  e_  ^'*1 

(45)  e  -e 

”^i  “2  “  f  (k) 


From  (43)  wo  have 
a^  (y-e^)  [Uq  (x,y)  -  u^  (>t,y-e^)  1 

(46) 

-a^<y)  [UQ(x,y+e^)  -  u^tjj.y)!  -  0 

If  we  take  aQ(x,y)  >  Uq(z)>  this  Is  trivially  true.  (¥e  must  Justify 
this  choice  in  subsequent  steps.)  And  (44)  simplifies  to 


(47)  eV.^'[a.  (y)  ^  (x)]  +  7^  'ta^(y)  (*«'/>!  “  0 


At  this  point  we  invoke  a  standard  device  in  homogenization  asymptotic 

analysis,  namely,  the  use  of  "correctors.”  We  assume 
d  3uq 

(48)  x^(y) 

k-l  K 

wlthX]^(*)  the  correctors.  Using  this  in  (47),  we  have  (agair 
using  the  summation  convention) 


(49)  V  Xk<»”  35^  *'S 3^*  ° 


If  we  take  X]^(y)  as  the  solution  of 

(50)  7^“(a^(y)  X,5(y)l  +  fa^^  <y-%)  -aj^(y)l-0 


(we  have  to  verify  the  well-posedness  of  (30)),  then  (49)  is 

f\j 

satisfied.  (The  term  u. (z)  is  determined  (formally)  from  the  0(e) 


term  in  the  system  (37)  (42)*) 


Regarding  the  eell-posedness  of  (50),  note  that 

(51)  7."[a.{y)  ♦(y)  1  -  ^  (y) 

d 

has  a  periodic  solution  on  Z  which  is  unique  up  to  an  additive 
constant  if  the  average  of  the  function 'i’  (y)  over  a  period  (eil)  is 
zero;  i.e., 

1  »• 

(52)  ?  ,  7(y+ke  )-  0  n-l,2;...,d 

t  )c*l  n 

This  condition  clearly  holds  izrt50),  and  so,  Xj^Cy)  is  well  defined 
(up  to  an  additive  constant). 


¥e  shall  determine  the  equation  for  Uq(z)  by  using  (48)  (50)  in 

(45)*  Using  the  Kroneeker  delta  function  5..,  we  have 
1  c- 

(53)  1  ^ 

-{  i  V  f  \  -V  ^*i  V  aTax, 

-  i  K 

-V '•i'lf'  x^(/)I  ♦  Uj]  -  f 

The  term  in  braces  is  zero  from  (50).  To  obtain  the  solvability 

condition  (52)  for  u^  in  (53),  we  introduce  the  average 


(54)  j  -  symmetric  part  (  {  -7^  [  (y)  Xj^  (f)  ]) 


Then  solvability  of  (53)  for  u^  gives  the  equation 


-».T  > 


(55) 


_1 

2 


d 

Z 

i»k"l 


’ii. 


f  (K) 


And  this  is  the  diffusion  equation  which  defines  the  limiting  behavior 
of  the  system  (36)  in  the  macroscopic  x-scale  in  the  limit  ase-l'O  . 


Ve  can  justify  the  asymptotic  cmalysis  by  using  energy  estimates 
or  probabilistic  methods  as  in  (Bensoussan  et  al.  1978).  (See  also 
Kunnemann  1983)*)  Ve  shall  omit  this  analysis  here. 


C.  Summary 


Returning  to  the  original  problem  (30)  for  the  evolution  of 
thermal  energy  on  a  microscopic  scale,  ve  have  shown  that  the  thermal 
density  u  (ttx)-*"  Uq  (t,x)  as  e+  0  (in  an  appropriate  norm)  where 


(56)  ^"0 

3t 


1 

2 


d 

Z 

i.j-1 


’ij 


^2 

3 


3x^3x. 


with 


(57)  **^3 


1 


I 


I 

Z  {V  "[a, 

k-1  ^  ^ 

+7j^"[aj^<y)Xj  (f)!) 


and  the  correctors  Xj^ ,  k  ■  1|2,...,d,  are  given  by 

Z  V“  la  (y)  X^  <y)  1  -  -  t  -®k  ^ 

(58)  i-1 


To  compute  the  limiting  "homogenized**  model  (36) »  one  must  solve  the 
system  (58)  (numerically)  and  then  evaluate  the  average  (57)> 
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Va  shall  assume  that  the  truss  has  a  regular  (e.g.,  triangular) 
cross-section  and  no  "interlacing*  supports.  ¥e  assume  that  the 
displacements  of  the  system  are  "small"  in  the  sense  that  no 
components  in  the  system  buckle.  ¥e  are  interested  in  describing  the 
dynamical  behavior  of  the  system  when  the  number  of  cells  (a  unit 
between  two  (triangular)  cross  sections)  is  large;  that  is,  in  the 
limit  as 

(59)  e:-  l/L  0 

¥e  shall  make  several  assumptions  to  simplify  the  analysis. 
First,  ve  shall  assume  that  the  triangular  sections  are  essentially 
rigid,  and  that  all  mobility  of  the  system  derives  from  the 
flexibility  of  the  members  connecting  the  triangular  ccmponents. 
Second,  we  shall  Ignore  damping  and  frictional  effects  in  the  system. 
Third,  we  shall  confine  attention  to  small  transverse  displacements 
Ti(t,x)  and  small  in  plane  rotations  ^(t,x)  as  indicated  in  Figure  9*2. 
¥e  shall  Ignore  longitudinal  and  out  of  plane  motions  and  torsional 
twisting.  Fourth,  we  shall  assume  that  the  mass  of  the  triangular 
cross  members  dominates  the  mass  of  the  Interconnecting  links. 

Systems  of  this  type  have  been  considered  in  several  papers 
including  (Noor  et  al.  1978)  (Nayfeh  and  Hefzy  1978)  (Anderson  1981) 
and  (Renton  1984).  In  those  papers  a  continuum  beam  model  was 
hypothesized  and  effective  values  for  the  continuum  system  parameters 
were  computed  by  averaging  the  associated  parameters  of  the  discrete 
system  over  appropriate  cell  volumes  or  areas.  Our  approach  to  the 
problem  is  based  on  homogenization-asymptotic  analysis  and  is  quite 


different  from  the  methods  used  in  these  papers. 

The  assumptions  made  above  simplify  the  problem  substantially. 
Qy  assuming  the  cross  sectional  components  to  be  rigid  and  ignoring 
out  of  plane,  longitudinal  and  torsional  motions,  ve  have  effectively 
eliminated  the  geometric  structure  of  the  truss.  ¥e  can  retain  this 
structure  by  writing  dynamical  equations  for  the  nodal  displacements 
of  the  truss  members.  For  triangular  cross  sections  nine  parameters 
describe  the  displacements  of  each  sectional  element.  The  analysis 
which  follows  may  be  carried  over  to  this  ease,  but  the  algebraic 
complexity  prevents  a  clear  presentation  of  the  main  ideas.  As 
suggested  in  (Noor  et  al.  1978)  one  would  need  a  symbolic 
manipulation  program  like  MACSYKA  to  carry  out  the  complete  details  of 
the  calculations.  We  shall  take  up  this  problem  on  another  occasion; 
for  now  we  shall  treat  the  highly  simplified  problem  which,  as  we 
shall  see,  leads  to  the  one  dimensional  Timoshenko  beam  (and  from 
there,  under  certain  constraints  on  the  parameters,  to  the  Rayleigh 
€uid  Buler  beam  models). 

I 

We  shall  begin  by  reformulating  the  system  in  terms  of  a  discrete 
element  model  as  suggested  in  (Crandal  et  al.  1980);  see  Figure  9.3. 

In  this  model  we  follow  the  displacement  rij(t)  and  rotation  4>^t) 

of  the  ith  mass  N.  The  bending  springs  )  tend  to  keep  the  system 

straight  by  keeping  the  masses  parallel  and  the  shearing  (k^)  tend  to 

s 

keep  the  masses  perpendicular  to  the  connecting  links.  We  assume 
small  displacements  and  rotations  so  the  approximations 


Flmre  5.3  A  lunped  parameter  model  of  the  sii^lified  truss  system. 


sin  (t)  Z  (t) 

(60) 

tan  (n^  (t)  A  :  (t)  /i 

are  valid. 


In  this  case  the  (approximate)  equations  of  motion  of  the  ith 


mass  are 


i  (t)  -  i  hi  {f (t)  j 

(61a)  ^  ^i+1 

t ^  1  > 


(The  spring  constants  depend  on  i  since  they  represent  the  restorative 
forces  of  flexed  bars,  bend  by  different  amounts.) 


.  (t) 

(61b)  n^Ct)  -  s"  {  -  ♦i(t)]  } 


where  we  have  normalized  a  ■  1  and  defined 
(62)  s"  s-  J  -  n^l 

and  similarly  for  S  (Ve  shall  omit  tz^ac  :er.  of  the  boundary 
conditions  at  the  ends  of  the  system.) 

To  proceed,  we  shall  Introduce  the  nondisMSional  variable  e  ■ 


and  rewrite  the  system  (61 )  as 


I1^(t)  -  -7^"  (Kg[V^%^(t)  -  } 


where 


,£+_  1 


v^'n.  -  7  (n.^^  -  n^),  1  <V'’i-i> 


Next  we  associate  a  position  in  the  system  Xj^€l[-1/2,  1/2] 

(normalizing  L  *  1)  with  each  mass;  and  we  Introduce  the  notation 


(65)  n(t,x^)  =rij^(t),  (|)(t,x^)  -  <}>^(t) 


Having  normalized  L  -  1 ,  we  have  c  ■  8.  and  x.  ,»x.  +  8.  -  x.  ♦£. 

i+1  1  1 

Let ■  {Xj^l  be  the  set  of  all  points  in  the  system.  In  this  notation 


(66) 


'\i 

xe  z 


(V  n)  (t,x)  ■-  [n(t,x+e)  -  Ti<t,x)l 
(V*””?!)  (t,x)  -  i  [n(tfX)  -  n(t,x-e)]  , 


and  the  system  is 

4.^(t,x,)  -K  (x,){7^V(t,x.)  - 

i  S  X  X  X 

(67)  +rV^'^{Kj^(x^)V^‘*^(t,x^)} 

Ti^(t,x.)  -  -V^‘’{k  (x  )  [7^‘'’Ti^(t,x. )  -  (t,x.)l} 

1  S  X  1  1 


The  scaling  of  (67)  may  be  interpreted  in  the  following  way: 

-2 

Formally,  at  least,  the  right  sides  of  both  terms  in  (67)  are  0( e  ). 
This  implies  that  the  time  variations  are  taking  place  in  the  "fast 
time  scale"  t  ■  t/e.  Also,  the  spatial  variations  are  taking  place  in 
the  microscopic  scale"  x  which  varies  in  e-increments  (e.g., 

*  t  )•  Introducing  the  macroscopic  scale  s  >  ex,  and  the  slow  time 
scale  o  "ex  ,  we  may  rescale  (67)  and  observe  its  dynamical  evolution 
on  the  large  space- time  scale  on  which  macroscopic  events  (e.g., 
"distributed  phenomena")  take  place. 


Rewritten  in  this  spatial  scale,  the  system  becomes 


-  .  i  K  (^)  f-) 

■>  e  s  e  e 


*i»  f 


(68a) 


dt 


-(()^(t,  ~  )}+ 


e 

1  ^e+r 


+  ^2 


(68b)  ^  -  i  6"-{k.  (^)  [6"V  (t.%]^ 


169;  6  “  eV  -  0(1)  in  e 


The  essential  mathematical  problem  is  to  analyze  the  solutions 


,  n^of  (68)  in  the  limit  as  e+0. 


B.  Mathematical  analysis 


To  proceed,  ve  shall  generalize  the  problem  (68)  slightly  by 

allowing  K  and  K  to  depend  on  z-as  well  as  z/c.  This  permits  the 
s  b 

restoring  forces  in  the  model  system  to  depend  on  the  large  scale 
shape  of  the  structure  as  well  as  on  local  deformations.  Ve  use  the 
method  of  multiple  scales;  that  is,  we  look  for  solutions  of  (tO)  in 
the  form 

n  (t)  ■  Ti^(t,z,y)  y=2/e 
<P  (t)  B  (p  (t,z,y) 


and  we  have 


(71)  ^  y  = 


Z  E— 

On  smooth  functions  ^(zf^  )  the  operators  ^  satisfy 
(6^'''  I/))  {z,y)  *  '*'  (z+E,y+l)  -'J' (Z/y) 

=  (z,y+l)  -  ili(z,y)  +  ip  (z+e,  y+1)  -  ip  (z,y+l) 

(72a)  +  3U)  1  3 

=  (ff  ip)  (z,y)  +E  (z,y+l)  +  -  - ^  (z,y+l)  +  0  (e'^) 

dz 


•  ••  •  •  -w*  ••  O  •  ^  • 


(z.y) 

(72\>) 


»  (z*y)  -  ’I'  (a-E»y“i} 

•  ip(z,y)  -'I'fe.y-i)  +’l'Cz,y-i)  (z-e,y-i) 


2 


(Z,y-1)  +0  (e  ) 


£  E 

We  assume  that  <1)  ,  n  may  be  represented  as 
ilJ^  (t,z,y)  =  ipQ  (t,z)  +  (t,z,y)  +  ... 

(73) 

(t,z,y)  =  tIq  (t,z)  +  £  rij^  (t,z,y)  +. . . 

and  substituting  (73)  in  (68)  and  using  (71)  (72),  ¥e  arrive  at  a 

sequence  of  equations  for  (  ^(j  ) » •  •  •  equating  the 

coefficients  of  like  powers  of  £ . 

-2-10 

Starting  with  £  ,  £  ,  £  ,  we  have 


(74)  i  s‘*'[r  K(z,y)  s"  \l)  (t,z)l  =  0 

which  is  trivially  true  from  (72b)  (73) •  The  same  term  for  (68b)  is 
trivially  satisfied  by  the  assumption  (73)*  Continuing 

(75)  £  ^  ^  \  ®  \  {s'*’nQ(t,z)  (t,z)  }] 

which  may  be  solved  by  using  the  corrector  X^(z,y)  and  taking 

(76)  (t,z,y)  =  X<j,  (z»y)  (t,z) 

with 

(77)  ®  (  (z,y)  s  Xa  (2,y))  ■  K  (z.y) 

If  we  regard  z  as  a  parameter  in  (77),  then  there  exists  a  solution 

98 


.  y.  •’v".  • 


>  .V 


,  unique  up  to  an  additive  constant,  if  K^(z,*)  are 

periodic  in  7,  if  there  exist  constants  A,B  so  that 

(78)  0  <A  £  Kjj  (z,y)  <  B  <  • 

and  if  the  average  of  K  (s,* )  is  zero 

S: 

L  ^  .i/2  *^8 

which  holds  if  the  system  is  pinned  at  the  ends  as  indicated  in  Figure 
5«3*  Let  us  assiue  that  (78)  (79)  hold,  and 

(80)  0  <  A  <  (z,y)  <  B  <  « 

(which  we  shall  need  shortly). 

Considering  (68b),  the  0( e~^)  term  in  the  asymptotic  expansion  is 

(81)  e  ^  ^  *a  ’*'0' 

Again  we  introduce  the  corrector  X^(z,y),  and  take  in  the  form 

(82)  rij^  (t,z,y)  -  '•'o 

which  gives  the  equation  for  the  corrector 

(83)  s'  {Kg(z,y)  [s'^x^  (z,y)  -Ij}  -0 
or 

(84)  s'{Kg(z,y)  “  K^(z,y)  -K^  (z,y-l) 

By  hypothesis  the  right  side  in  (84)  is  periodic  in  y  and  has  zero 


average  (79}>  Hence,  (84)  Has  a  periodic  solution,  unique  to  an 
additive  constant. 


Continuing,  the  0(c)  term  in  (68a)  is 
s'^  {r  Kjjte.y)  s>2 
+Kg  (z,y)  ts'*'  (t,2,y)  ip^(t,z,y)] 

3  0 

+K  (z,y)  -5—  (t,z) 

S  dZ 

+s*  {r  (z,y)  (t,2,y)} 

+S'*'  {r  Kjj  te.y)  —2  ’^0 
dz 

+  —  {  r  K.  (Z,y+1)}  —  'I'q  (t**) 

3z  3*  2 

2  ^  ’^C 

+  — 2  {  r  l^(z,y+l)l 
3z^ 


This  should  he  regarded  as  an  equation  for  as  a  function  of  7  with 
(t,z)  as  parameters.  In  this  sense  the  solvability  condition  is  as 
before,  the  average  of  the  sum  of  all  terms  on  the  left  in  (83) > 
except  the  first,  should  be  zero.  ¥e  must  choose  so  that  this  in 
fact  occurs;  and  that  defines  the  limiting  system. 


Using  the  correctors  (76)  (82),  we  must  have 


Average 


3^^ 

—5^  [  s'^(r  (Z,y)+  s'^  (r  Kjj  (S,y)  te.y)  1 

3z^ 


(e6)  - 

3z 

+  Kg  (z,y)  is*  fe»y)  -  X^fc»y))>  -  0 
Defining  the  functions  Bl(z),  6(z)  by  the  associated  averages  in  (86), 


the  averaged  equation  is 


¥e  have  shown  that  a  siaplified  aodel  of  the  dynamics  of  the 
truss  with  rigid  cross  sectional  area  may  be  well  approximated  by  the 
Timoshenko  beam  model  in  the  limit  as  the  number  of  cells  ('^LA) 
becomes  large.  The  continuum  beam  model  emerges  naturally  in  the 
analysis,  as  a  consequence  of  the  periodicity  and  the  scaling. 

To  compute  the  approximate  continuum  model,  one  must  solve  (77) 
and  (84)  (numerically)  for  the  correctors  and  then  compute  the 
par^uDeters  in  (87)  (88)  by  numerically  averaging  the  quantities  in 
(86)  (and  its  analog  for  (68b))  which  involve  the  correctors  and  the 
data  of  the  problem. 


6.  HoBog«nisation  and  Optimal  Stochastic  Control 

In  this  section  and  the  next  we  show  that  the  process  of  deriving 
effective  "continuum"  approximations  to  complex  systems  may  be 
developed  in  the  context  of  optimal  control  and  state  estimation 
designs  for  those  systems.  This  procedure  is  more  effective  than  the 
procedure  of  first  deriving  homogeneous/continuum  approximations  for 
the  structure,  designing  a  control  or  signal  processing  algorithm  for 
the  idealised  model,  and  then  adapting  the  algorithm  to  the  physical 
model.  In  fact,  separation  of  optimisation  and  asymptotic  analysis 
can  lead  to  incorrect  algorithms  or  approximations,  particularly  in 
control  problems  where  nonlinear  analysis  (e.g.,  of  the  Bellman 
dynamic  programming  equation)  is  required.  The  problems  treated  here 
and  in  the  following  section  are  abstract  systems  which  illustrate  the 
basic  techniques.  At  the  beginning  of  section  7  we  shall  present  a 
simple  argument  which  shows  how  the  class  of  models  treated  here  might 
arise.  In  subsequent  work  we  shall  apply  the  combined  homogenization 
-  optimization  procedure  described  here  to  the  problem  of  controlling 
the  dynamics  of  lattice  structures  like  the  truss  structure  analyzed 
in  the  previous  section. 

6.1  A  Prototype  Problem 

The  interaction  of  homogenization  and  stochastic  control  was 
discussed  briefly  in  the  book  (Bensoussan,  Lions  and  Paponicolaou 
1978),  and  in  (Bensousscui  1979)  and  (Blankenship  1979).  The  recent 
paper  (Bensoussan,  Boccardo  and  Murat  1984)*  provides  the  first 
systematic  analysis  of  an  abstract  control  problem  involving 
homogenization.  Ve  shall  briefly  summarize  its  main  results  against 


the  background  of  the  lattice  syatwB  dlacuased  in  section  5 


The  problem  is  to  control  z  (t)  given  by 

dx^  -  [  g  J  i  b  I  x^)]  dt  +  o  tx^,  j  x^)  dw  (t) 

(1) 

x^  (0)  -  X,  0  <  t 

with  z^(t}  defined  in  a  bounded  domain  OCR  with  smooth  boundary. 
Here  g,  b,  and  a  are  smooth  functions  of  their  arguments,  w(t)  is  a 
standard  R°  -  valued  Wiener  process,  v  is  the  control  and  e>  0  is  a 
parameter.  We  assume  that  g,  b,  and  a  are  periodic  in  their  second 
argument  with  period  one  on  the  unit  torus  Y  in  R  . 

Let  T  be  the  first  exit  time  of  z  (t)  from  the  domain  0.  The 
coat  function  is 


(2)^  (vO) 


■“Jo' 


L(X®, 


i 


exp 


X®,  ^  ,v)ds)  dt] 


and  we  define 

u^(x)  -  inf  (v(*)) 
v(*) 

We  assume  that  the  cost  rate  L(z,y,v)  is  periodic  in  y  on  the  torus, 
and  has  linear  growth  in  v.  The  discount  factor  c(z,y,v)  is  uniformly 
bounded,  positive,  cmd  periodic  in  y.  The  set  of  admissible  controls 
U  consists  of  feedback  functions  v( • ) 

(4)  V  (t)  -  <|)^  (x^  (t) ,  j  (t)  ,t) 


*  We  are  grateful  to  Professor  A.  Bensoussan  for  transmitting  a 
preprint  of  this  paper  to  us. 


(¥e  shall  Justify  this  presumed  structure  for  the  control  law  in  more 


detail  later.)  The  system  has  highly  oscillatory  coefficients;  and, 
as  in  section  3>  one  would  expect  the  solutions  of  the  control  problem 
to  be  well'  approximated  by  the  corresponding  control  piroblem  with 
g(x,y,v),  b(x,y),  a(x,y),  L(x,y,v),  and  c(x,y,v)  replaced  by  their 
averages  (appropriately  defined)  over  y.  This  is  not  precisely  the 
case  and  one  must  carry  through  a  complete  asymptotic  analysis  to 
determine  the  exact  nature  of  the  limit  and  the  approximation. 

In  (Bensoussan,  Boccardo,  and  Murat  1984)  this  analysis  was 
carried  out  in  terms  of  the  Hamllton-Jacobi-Bellman  (HJB)  equation  for 
the  optimal  cost,* 

A^u  ^  =  H  (x,  i  X,  u^,  Du^) ,  u^lp  *  0 


(6)  »  -  a.  .  (x,— ) 


X,  a  1  ^  ,  x,  3 

ij  '"'e'  3x,3x,  ”  e  ‘’i  3x. 

i  j  1 


where  a  .  .  ■  i  (  ao^  )  .  .  emd 
iD  2  '  ' 


(<7)  x,y,q,p)  =  inf  {l  (x,y,v)  +  p*g  (x,y,v)  "  q  c  (x,y,v)} 


Notice  that  the  Hamiltonian  H  is  periodic  in  y.  The  objective  of  the 
analysis  is  to  determine  the  limit 

u(x)  =  lim  u^(x) 

(8)  e  0 


*  Here  and  in  the  following  we  use  the  summation  convention  that 
repeated  indices  iJ  are  summed  over  their  full  range. 


to  identify  u(z)  as  the  solution  to  a  HJB  equation,  and  through  this 
to  associate  a  "liniting  stochastic  optimal  control  problem"  with  the 
original  problem.  The  method  is  "homogenization”  of  the  nonlinear 
partial  differential  equation  (5).  To  accomplish  this,  it  is 
necessary  to  impose  further  regularity  and  growth  conditions  on  the 
coefficients  in  (5 )-(?)• 

Assume  U  .  is  a  nonempty  subset  of  a  compact  metric  space  U;  and 
ad 

(9)  L(x,y,v)  j  r"  X  r”  X  -►R 

is  continuous,  periodic  in  y,  and 

(10)  m^lvl^  -  m^  ^  L  (x,y,v)  ^  h  (  l+lv|^  ),  ^  0 

Also, 

9(x,y,v):  r\  r"  x  U 

(11) 

a(x,y,v)!  r"x  r"  x  U  -►r"*" 

ad 

are  periodic  in  y,  continuous  in  their  other  arguments,  and  satisfy 

|g| •<  i  (  i+lv|) 

(12) 

0  <  c  c 

Under  these  assumptions  standard  selection  theorems  guarantee  the 
existence  of  a  control  law  v(z,y,q,p}  achieving  the  inflmum  in  (7)  and 
for  any  q,p  fixed 


I  » 


(13)  ^  ^  (v,y,q,p)  I  <  c(  1+  1p1  +  ) 

I  H(x,y,q,p)|  £  c(  1+  |p|  ^  +  kl) 


Moreover, 


H{  x,y,q,0)  <  c  (  1+  Iv]^)-  q6»  <1^0 
H(x,y,q,0)  >  -Cj^  +  qc-q6,  q  £  0 


for  some  non-negative  constants  c,  and  some  arbitrary  v  e  ^ad 
Also , 


i 


i  t 


I  C 


B(x,y,q,p)-  H(»,y,s,t)  <  (ft)  •  t  (x,y,0(x,y,q,p)) 

(15)  -(q-s)  c (x,y,v(x,y,q,p) ) 

<  c|  p-t  I  {1+  |p|  +>|  kh-tcla-sl 

and  a  similar  condition  holds  with  the  roles  of  (q,p)  and  (syt) 
reversed.  These  growth  conditions  suffice  for  the  asymptotic  analysis 
of  (5). 


6.2  Invariant  Measures  and  Correctors 


Consider  the  second  order  operator 


:i6)  A  =  A^-  -a^j(x,y) 


-  b.  (x,y)  3— 


and  its  formal  adjoint 


(17)  A  -  -  I7-  <  *ij  ItT^  ^  ^ 


For  any  x  fixed  let  aCx.j)  be  the  solution  of 

* 

A  iw-0,  m(x,y)  periodic  in  y 

(18) 

m  >  0,  f  m(x,y)dy  -1 
Y 

2  o 

ifith  m  regular  (in  ,  2  ^  p  <■  *).  Since  x  is  restricted  to  a 

compact  subset,  we  may  assume 

(19)  ®  ^  2  £  ® 

for  some  constants  m  and  m.  Thus,  m(x,y)  defines  a  probability 
measure  on  Y  for  each  x,  which  we  call  the  invariant  measure 
associated  with  A. 

A  key  assumption  in  the  method  is  that  the  drift  term  b(x,y)  in 
(l )  .is  "centered*  in  the  sense 

(20)  /y  m(x,y)  b(x,y)  dy  »0 

If  this  assumption  fails,  then  the  asymptotic  analysis  takes  a  very 
different  form  from  what  follows,  and  the  results  (i>e.,  the  limits) 
have  a  totally  differ,ent  character. 

The  centering  hypothesis  and  the  regularity  assumptions  mean  that 
the  "correctors"  defined  by 

\  (x*y)  *  (x,y)  ,  i-1,2,...,  n 

(21)  . 

y  X  (Xfy)  periodic,  x^  (x,y)  dy  -0 
exist  and  are  smooth  (l.e.,  they  are  functions).  These  functions 


play  a  key  role  in  the  homogenization  procedure  (Bensoussan,  Lions  and 
Papanicolaou  1978)  in  defining  the  limiting  system.  Solution  of  the 
system  (21 )  is  the  key  numerical  problem  in  applying  the 
homogenization  techxiique. 

6*3  Identification  and  Interpretation  of  the  Limit  Problem 
Using  the  invariant  measure  and  the  correctors,  we  define 


Q..  (x)  «  /  m(x,y)  [a  (x,y)  -  a.. 
^  y  ij  ik 


(22) 


“  s . ,  3y  1  T  i 

- 2  x^+bjX  )]<3y 


/y  m(x,y)  [  a.  .  (x,y)  -a^,^  ^  ]dy 

(23)  A  u  -  -  -  r .  (X) 


^24)  H  (x,q,p)  =/  m(x,y)  H{x,y,q  (I-Dx)p)  dy, 

^  .  k 

(DX  )  =  £X_ 


P  3 


Using  the  definition  (21)  of  the  correctors,  we  can  rewrite 


Q^^(x)  as 


(25) 


Qij  (X)  -  m(x.y)  [a.^ 


j 


-  a . 


jk  9y^ 


+a 


rk  3y 


ixL  , 


>v. 


> 
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which  is  uniformly  positive  definite*  This  with  the  other  assumptions 
means  that 


,  ^  A  u=  H  (x,u>Du)  ,  u L  -0 

(26)  -  -  r  ^ 

has  a  unique  solution  in  W  ,  pe  [2,  ®) 


C: 


Equation  (26)  defines  the  "limiting  control  problem"  associated 
with  (1)  -  (4). 

THEOREM  (Bensoussan,  Boccardo  and  Murat  1984)  Under  the  assumptions  of 
regularity  and  nondegeneracy  of  a^^ (x,y) 


i  I 


(27)  -►  u  weakly  in 


for  some  pQ,>  2« 


i  c; 


The  function  H  in  (24)  may  be  rewritten  as 


fz 

i'' 

*  S 

t  •*  »  ^ 


•  •  *  * 

I  •  ' 

t  t: 


H  (x,q,p)  *  inf  {  /  in(x,y)  [L  (x,y,v(y)) 
+  (g^  (x,y,v(y))  -  3^  gjj  (x,y,v(y)) 


(28)  -<J  c<x»y,v(/))]  dy| 


A  inf  {  L  (x,v(*))  +  P  •  g  (x,v  {•)) 

V  (•)  , 

-q  c  (x,v  (•))/ 


u*  • 


t  c 


From  the 


•I’  *A •  j-  -*  .* . 


.'I*.'-'.  If.  -T 


where  v(*)  is  any  Borel  function  on  Y  irith  values  in 
final  expression  it  is  clear  that  H  is  the  Hamiltonian  of  a  control 
problem. 

Using  (22)  (23)  and  (28),  we  can  identify  this  problem  explicitly 
as 


(29a)  'i(x)  -  inf  J  (v  (•)) 

v(*) 


(29b) 


{v{*))  ■  E^(.)  C/q  L  v(t)) 


exp  (-/q  c  (x, (s),  v(s))  ds) 


dt] 


dxa 


[g 


(29c) 


(x  (t) ,  V  (t)  +r  (x  (t) )  ]  dt 

+  '^2Q  (X  (t) )  dw  (t) 
X  (0)=x,  0  <  t 


with  the  first  exit  time  of  the  (controlled)  process  x(t)  from  the 
compact  domain  0. 


A  key  property  of  the  limiting  control  problem  is  that  the 
admissible  control  laws  depend  on  y,  the  "rapidly  varying  state 
variable",  as  may  be  seen  from  (28).  Thus,  the  "fine  structure"  of 
the  original  problem  (1)  -  (4),  that  is,  the  periodic  dependence  of 
controls  on  y^(t)  ■  x^(t)/£  ,  is  retained  in  the  limiting  problem.  In 
effect,  the  limiting  optimal  control  law  depends  on  the  fine  structure 
-  property  which  may  not  be  desireable  in  some  engineering 
implementations . 


no 


Notice  that  the  limiting  state  dynamics  (29c)  emerge  naturally 
from  the  asymptotic  analysis  (see  Bensoussan,  Boccardo  and  Murat 
1984  for  details)  of  the  nonlinear  HJB  system  (5)  -  (7)»  Note 

further  that  simply  averaging  the  functions  g,  h,  a  ,  L,  c  with 
respect  to  y,  and  then  posing  an  optimal  control  problem  in  terms  of 
the  averaged  (y  Independent)  functions  leads  to  wrong  answers  for  two 
reasons.  First,  the  appropriate  averaging  process  involves  the 
invariant  measure  and  the  correctors,  and  the  role  of  the  latter  is 
not  obvious  in  a  naive  application  of  averaging  methods.  Second,  as 
(29)  shows,  the  optimal  control  law  that  emerges  in  the  limiting 
process  depends  on  y,  which  cannot  be  the  case  when  the  averaging 
process  is  separated  from  the  optimization  process. 

The  key  numerical  problems  in  applying  the  homogenization  -> 
optimization  procedure  to  a  specific  problem,  e.g.,  control  of  the 
lattice  structure  described  in  section  are 

1.  solving  for  the  invariant  measure  (18), 

2.  verifying  the  centering  hypothesis  (20),  solving  for  the 
correctors  (21), 

3.  computing  the  averaged  quantities  (22)  -  (24),  and 

4.  solving  the  limiting  control  problem  (26)  -  (29). 

Sequential  solutions  of  these  problems  constitute  algorithms  for 
simultaneous  homogenization  and  control. 


6*4  Application:  Homogenization  -  optimization 
of  lattice  structures 


To  illustrate  the  techniques  of  the  last  sections,  we  shall 
reconsider  the  model  for  the  lattice  structure  analyzed  in  subsection 
with  control  actuators  added.  The  truss  shown  in  Figure  3.3  is 
again  constrained  to  move  in  the  plane  and  torsional  motion  is 
excluded  to  simplify  the  model  and  confine  attention  to  the  basic 
ideas.  Now,  however,  we  include  a  finite  number  of  actuators  acting 
to  cause  transverse  motions.  !^e  truss  with  actuator  forces  indicated 

9 

by  arrows  is  shown  in  Figure  6.1.  The  corresponding  discrete  element 
model  is  shown  in  Figure  6.2. 


Suppose  that  the  physical  actuators  act  along  the  local  normal  to 
the  truss  midline  as  shown  in  the  figures,  and  that  the  forces  are 
small  so  that  linear  approximations  to  transcendental  functions  (e.^., 

f\j 

sin  etc.)  are  valid.  Then  the  controlled  equations  of  motion 

of  the  discrete  element  system  are  (recall  equation  (3.63)) 


r  ^  (t)  -  i  ^  t  ^  (t)  ]  +V^'^[  kJ  (t)  ] 


n  ^  (t)  -  -V^"  {  (t)  -  (t)]}  +  E  5  Ci,i.)  u. 

3  X  i  3  D 


where  the  notation  in  (5*64)  has  been  used. 


Figure  6.2  Discrete  element  model  of  the  controlled  truss 


Hence ,  if 


and  are  the  locations  of  the  actuators. 

6(ifi  J  *  0  for  all  there  is  no  actuator  located  at  the  ith 

point  which  correponds  to  the  physical  point  e  [0,l].  The  number  m 
of  actuators  is  given  at  the  outset  and  does  not,  of  course,  vary  with 
the  scaling. 

The  control  problem  is  to  select  the  actuator  forces  as  functions 
of  the  displacements  and  velocities  of  components  of  the  structure  to 
damp  out  motions  of  the  structure.  Measurements  would  typically  be 
available  from  a  finite  number  of  sensors  located  along  the  structure. 
Ve  shall  not  elaborate  on  this  component  of  the  model,  and  shall 
Instead  assume  that  the  entire  state  can  be  measured.  To  achieve  the 
stabilization,  we  shall  associate  a  cost  functional  with  the  system 
(30).  Let 

(32)  “  [u^  (t) ,  . . . ,  (t)  1 

be  the  vector  of  control  forces,  and 

[  u  01  -  /q  [ml  (t)]^  +  b.  [  (t)]^  +  oi.[  ;I»  ^ 

in 

+3^  [  (t)]  +  Z  5  (iyi.)  u?  (t)}  e  ^  dt 

j-i  3  3 

where  (a^.,b^  )  and  ( a^,  6^)  are  non-negative  welgnts.  Formally,  the 
control  problem  is  to  select  6(i,ij)u^(t),  i-1,...,N,  j-1,...,m  to 
achieve 


(34)  infJ^IuOl 


subject  to  (30)  (31)  and  the  appropriate  boundary  conditions.  The 
case  0  corresponds  to  stabilization  by  feedback. 

The  analysis  of  this  control  problem  is  based  on  the  scaling  used 
in  section  5»  equations  (5«63)  -  (5»69).  Let  t  ■  t/e  be  the  fast 
time  scale,  then 


V  ^  ^  2  '\j  £  ^  2 

lu(*)]  ■  /q  £  1  {a^[  ^  (T)  ]  +b^  [  n  ^  (T)] 

i=l 


+  a  £^  [  ij;  ^  (T)^  +e  £^[  n  J  (t)]^  +  E  6  a,i.)  [u.  (f)  ]^}e"''^‘  dT 

j-1  J  J 

with  (t)"  (£T)/etc. 

Let  (  'J' ,  'P  ,  n  ,n  )  be  the  state  vector  of  the  system  (30)  with 
'J'  •  “d  similarly  for  the  other  terms.  Let 

V  •  •  'P',  T1  ,Ti  )  be  the  optimal  value  function  for  the  problem 

(30)  (35).  Then  the  Bellman  problem  associated  with  (30)  (35)  ia 


i2i  -ere 


^  ii '*1 '♦i  *  v 

.  e  !  (  t  t  n,  :  4.,!  *  i  ♦J)  V  ♦, 


i=l  r 


i»l 

min  {  £  I  2  [  «j  ^  ^ 

u.  j  i-1  ->  ■»  i 


.e  ?  [  a^Tj.b^n^e^  -eyv.o 

i*l 


VZw- -W.v 


BENAfiKS: 


1.  Note  that  the  minimiaation  in  (36)  is  well  defined  if  the  admissible 
range  of  the  control  forces  is  convex  since  the  performance  measure  has 
been  assumed  to  be  quadratic  in  the  control  variables  d(ifi^)u^  • 

2*  Since  we  have  not  included  the  effects  of  noise  in  the  model »  the  state 
equations  are  deterministic  and  the  Bellman  equation  (36)  is  a  first 
order  system.  To  "regularize"  the  analysis,  at  least  along  the  lines 
followed  in  conventional  homogenization  analysis,  it  is  useful  to 
include  the  ef'fects  of  noise  in  the  model  and  exploit  the  resulting 
coerclvity  proper! tes  in  the  asymptotic  analysis. 

3.  If  we  introduce  the  macroscopic  spatial  scale  z  ■  ex,  the  mesh  {x^  }, 
and  the  variables 

(37)  ip  (t,z^)  «  (t),  ^  (t,z^)  -  'i'  ^  (t),  etc. 

then  the  sums  e  I  msty  be  regarded  as  Rlemann  approximations  to 
i 

Integrals  over  the  macroscopic  spatial  scale  z.  The  asymptotic 
analysis  of  (36)  with  this  interpretation  defines  the  mathematical 
problem  constituting  simultaneous  homogenization  -  optimization  for 
this  case. 


¥e  shall  return  to  this  challenging  problem  in  subsequent  work. 


7*  HofflogenlzatioQ  and  State  Estimation  in 
Heterogeneous  Structures 

7.1  Problem  Statement  and  Background 

Signal  processing  problems  arise  in  the  control  of  large  space 
structures  in  several  ways.  Our  special  Interest  here  is  in  the 
treatment  of  a  nonlinear  filtering  problem  for  a  prototype  abstract 
jystes  with  a  homogeneous  infrastructure.  ¥e  shall  give  ft  detailed 
treatment  of  the  filtering  problem  fqr  the  system 

dx^(t)  -  g(^^-^ldt  +  o[^^^^ldw(t) 

dz'^Ct)  -  h[^^^^]x*^(t)dt  +  dv(t) 

x^O)  -  C  .  z^O)  -  0.  0  <  t  <  T  ,  E  >  0 

where  C  ia  an  r"  -  valued  random  variable,  g,  a  ,  and  h  are  periodic  on 
the  (unit)  torus  in  R  ,  and  w(t)  and  v(t)  are  Independent,  standard 
vector-valued  Wiener  processes  which  are  independent  of  ^  .  The 
filtering  problem  for  (l)  is  to  estimate  x^(t),  i.e.,  compute  its 
conditional  density,  given  “a|z^(s),  0^  s_£  t},  the  o-algebra  of 
observations.  We  are  interested  in  the  behavior  of  this  filtering 
problem  in  the  limit  as  £  0. 

In  the  model  (l)  the  vector  x^(t)  may  be  regarded  as  the 
composite  state  of  the  overall  system  formed  from  the  lexicographical 

*  The  results  in  this  section  are  Joint  work  with  A.  Bensoussan 
at  INRIA  in  Versailles.  This  research  was  also  supported  in  part  by 
the  Department  of  Energy. 


listing  of  the  states  of  each  of  the  components  of  the  sgrstem.  The 
periodicity  of  g,  a  ,  and  h  represents  a  regularity  property  of  the 
array;  the  small  parameter  e  represents  a  natural,  non-dimensional 
"distance"  or  "coupling"  variable  characterizing  component 
interactions.  In  a  subsequent  paragraph  we  shall  describe  a  prototype 
system  in  this  class. 

One  would  expect  the  system  (l)  to  be  well  approximated  as  £  0 
by  a  similar  system  with  g(x/£),  a(x/e),  and  h(x/£)  replaced  by  their 
averages  g,  T ,  and  h  over  the  torus.  This  is  the  case,  although  the 
precise  nature  of  the  average  is  difficult  to  guess  from  a  cursory 
inspection  of  ( l) .  The  filtering  problem  for  the  associated  limiting 
system  is  just  the  Kalman-Bucy  filtering  problem  which  has  a  simple, 
closed  fozm  solution.  By  constructing  an  asymptotic  expansion  for  the 
conditional  density  of  x  (t)  given  ,  we  can  obtain  a  family  of 
finite  dimensional  linear  filters  which  (presumably)  provide 
increasingly  accurate,  e.g.,  0(e),  0(e  ),...,etc.,  approximations  of 
the  conditional  density  of  x^(t).  The  technique  used  to  derive  the 
result  is  "homogenization"  of  a  linear  stochastic  partial  differential 
equation  for  the  (unnormalized)  conditional  density  of  x^(t)  given  Z  . 

While  the  system  (1 )  is  obviously  only  an  example  of  a  larger 
class  of  problems,  we  shall  see  that  its  analysis  has  all  the 
essential  difficulties  of  more  general  problems.  Before  starting  the 
analysis  it  is  useful  to  illustrate  how  a  problem  like  (l)  might 


Consider  the  prototype  system: 


1  N 

dx,(t)  •  a[Xj(t),  u(t)]dt  I  b[x.(t)]dw  (t)i  “  N, 

i  i  J 


a  and  h  are  smooth  functions  of  their  (vector- valued)  arguments, 

and  are  standard  (vector)  ¥iener  processes  which  are 

lendent  for  (i,j)  *  (k,®-)  and  u(t)  is  a  vector  of  control 

iblea.  The  functions  a  and  b  are  the  same  for  all  the  subsystems 

T 

the  overall  system  with  state  x(t)  ■  [x2-(t),...,XQ(t)]  has  a 

;eneous  structure.  The  coupling  is  random  and  normalized  by  l/N 
fleet  the  assumption  that  each  subsystem  has  0(l)  coupling  to  the 
nder  of  the  system  (as  opposed  to  0(n),  0(1/N),  etc.),  no  matter 
arge  the  latter  is. 


Associated  with  (2),  we  define 

N. 

S(t)  »  E  X. (t)  “  "the  aggregate  output" 
1=1  ^ 

1  ” 

o(t)  "  «  E  X, (t)  =  "the  average  output" 

1=1 


>se  that  in  the  process  of  controlling  the  system,  we  observe  not 
.  but  the  aggregate  S(t)  through  the  measurement 


dz(t)  =  h[S(t)ldt  +  dv(t) 


h  smooth  €Uid  v(t)  a  standard  Wiener  process.  Suppose  further 
the  control  u(t)  is  defined  by  u(t)  =  f[s(t)]  with  S(t)  an 
mate  of  S(t)  derived  from  z(s),  s  <  t.  We  would  like  to  analyze 
(4)  in  the  limit  as  N"*'*  ;  and,  more  precisely,  to  show  that  this 


analysis  Involves  the  asymptotic  analysis  of  systems  scaled  like  (l) 


N 


Defining  (t)  ■  x^(t)  -  a(t),  we  have  “  0*  The 


aggregate  output  S(t}  satisfies 
N 

dS  (t)  »  I 


(5) 


i-1 


a  (>c^  (t) ,  u  (t)  )dt 


1  ? 
j-i 


M 


b  {x .  (t) )  S 


i=l 


dw  .  (t) 
iD 


=Na  {p  it)  ,u  (t)  )dt 


2  1  ** 

+0  (5 1 X .  (t)  1  )  dt+b  to  (t)  2  (t) 

^  i,j=l  ^ 


N 


N 


+b  to  (t)  z  6x .  (t)  Z  dw  .  (t) 
*  **  i-1  ^  j=l 


+0(|6x^{t)|^)dw(t) 


where  w(t)  is  a  vector  Vlener  process  defined  from  the  components  of 


2 

w^j(t).  Neglecting  0(1 6 x^(t) 1  )  terms,  we  have 


(6) 


do(t).  =  a(a(t),  u(t))dt  +  b(<T(t))— »  .Z  dw  (t) 

i,j.l 


To  treat  the  last  term,  we  use  the  formal  argment  in  (Seman  1982) 

which  goes  as  follows:  As  N'*’*  a  "local  chaos"  condition  prevails  in 

which  each  subsystem  with  state  6x^(t)  behaves  "independently"  of 

N 

every  other  subsystem,  and,  in  effect,  of  the  noises  (S  dWx^(t)/H), 

i-1 

1~1,...,N.  That  is,  a  law  of  large  numbers  applies  to  the  last  term 


as  N  -*«> .  Since 


N 

I  6x.  (t)-0 


W.  ■.  ■»  .•*  A  TV  .■»  .-B  Tf 


r: 


I  I 


-.  .V 
;■  ■^• 


i 


z 


i  t 


I  r; 


■*  I'. 

»  L 


t  c 


by  the  definition  of  a(t),  the  last  term  in  (6)  is  sero.  (In  a  more 
general  situation,  this  term  would  approach  sero  as  N  •»■")  Notice  that 


1  ** 

-  I 

i.j=l 


w.  .  (t) 


in  the  second  term  behaves  like  a  standard  Wiener  process  for  each 
Thus,  for  large  N  we  obtain  the  approximate  model 


(t) 

do(t)  -  a(o(t),u(t))dt  +  b(o(t))dw(t) 


Now  let  a(a,S)  ■  a(a,S)  and  assume  that  S  and  S  have  the  same 
order  behavior  in  N  for  N  large.  Defining  e  ■  l/N,  we  have  two 
descriptions  of  the  aggregate  behavior  of  (2)  for  N  large 


(8a) 


(8b) 


do(t)  -  a(o(t),  ■^o(t))dt  +  b(o(t))dw(t) 

dx(t)  -  h(  ”0  (t))dt  +  dv(t) 


dS(t)  .  a(cS(t),S(t))dt  +  i  b(eS(t))dw(t) 


dz(t)  *  h(S(t))dt  +  dv(t) 


So  to  analyze  the  aggregate  behavior  of  the  original  system  (2)  as 
N-*<»  ,  we  can  study  (8a)  or  (8b)  as  0.  If  a,  b,  and  h  have  a 
periodic  or  randomly  recurrent  dependence  on  their  arguments,  then  the 
analysis  of  (8a, b)  involves  a  homogenization  problem. 
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The  literature  in  mathematical  physics  and  engineering  contains 
many  examples  of  systems  scaled  like  (8)  which  can  be  effectively 
treated  using  homogenization  theory  (Keller  1977)  (Larsen  1975t  1976). 
Homogenization  methods  have  not  been  developed  in  control  theory, 
other  than  the  brief  treatments  in  (Blankenship  1979)  (Bensoussan 
1979). 

7.2  Preliminary  Analysis 

Let  (  ^ ,P,P)  be  a  probability  space  on  which  are  defined  two 
Independent  Viener  processes  w(t)  and  z(t)  with  values  in  and 
respectively.  Let  C  be  a  Gaussian  random  variable  with  values  in  r" 

'Xj,  . 

which  has  mean  x^  eind  covariance  Pq.  Suppose  ^ is  independent  of  w(t; 

4a  0t> 

and  z(t).  Let  P  ,  t  ^  0,  be  a  family  of  o-algebras  with  P  ■  P,  such 
that  w(t)  <md  z(t)  are  adapted  to  and  C d.s  -  measureable.  Let 
Z  *'a{z(s),  s<t}.  Let  Y  be  the  unit  torus  in  R*^  and 

g(y)  c  L  (r";  r") 

(9)  o(y)  E  L  (r”;  R*')  ;  invertible 

h(y)  e  L  (r";  R^) 

which  are  defined  on  the  torus  Y,  and  which  are  sufficiently  smooth 
there . 

Let  x^(t)  be  the  solution  of  the  I to  equation 


and  note  that  x^^Ct)  is  independent  of  z(t).  Consider  the  processes 


c 

.*/ 

I 


(11) 


w^(t)  “  -  a  g(— )x^ds  +  w(t) 

0  ^ 


v^(t)  “  -  h(— )x^  ds  +  z(t) 

-  e 


and 


(12) 


€  ^  e 

U^(t)  “  exp{  h(— )x^*dz  +  o  g  (— )x^’d\^ 

0  .  0  ^ 


h  Ih  (^)x^l?  ds  -  *s  |a  ^g(^)x^l^ds 
0  e  u  e 


For  any  finite  T,  one  has 

(13)  E  m®^(T)  <  « 

which  is  a  consequence  of  the  following  condition  (see  A.  Bensoussan, 
J.  L.  Lions  1978) 

(14)  E  exp  6|x^  (t)  1^  ^  c,  V  t  e[o,T] 


1  _  * 


To  check  (14),  consider  the  backward  Cauchy  problem  (a  “  "J  , 
and  we  shall  use  the  summation  convention  from  here  on) 
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»  J 


r:  ■  •  — 

I 


V»,' 


+  a  ^  ^ 

-  0  ,  s  ^  t 

(15) 

as  Ij  e  oX^oXj 

u(x,t}  ■  exp 

(6|x|^) 

Then 

(16) 

E  exp  ( 6 |x^(t) 1^  )  - 

E  u(C.O). 

Consider  the  function 

(17) 

t(x,8)  =  exp[P(s) 1 

x]^  +  pfe)  ], 

P(s)  _>  0 

P(t)  =  6  . 

P(t)  -  0 

¥e  have 

• 

3^  +  a 

33  i j  3x. 3x . 

.  .  ^  -jJ 

■  =  a  P  +  p|x|^ 

*y  0 

+  2tr  a  P+4lax|  P^] 

(18) 

<  a  P  +  (P+4|  |a| 

+2ltr  a|p] 

<  C(P  +  C'+4l  |a|  jp^) 

+2n  1 laj |pl 

Choosing  P  and  p  so  that 

(19) 

P  +  4llallP^  =  0, 

p  +  2n  1  lal  Ip  =» 

0 

we  have 

(20) 

P  (s)  =  6/[  1-4|  |al  |6 

(t-s)  ] 

expfe)  =  l/[  l-4[  |al  [5 

By  the  maximum  principle, i 

%  (x,a)  >u(x,s). 

Hence , 

E  exp  «lx^(t)|  <  E  [ 

(1  -  4  Hall  6t)"'2 

(21 )  1 1  1 1  —12 

exp6[(l-4l|  all  6t)I  "  2Pq  ]  ^jx^l 
/|(l-4l|a||6t)I  -  26Pq| 

Therefore,  sufficient  conditions  for  (14)  to  hold  are 


1  -  4  II  a  II  6T  >  0 

(22) 

(1-4  II  a  II  6T)  I  >  26Pq. 

which  hold  if  5  is  sufficiently  small.  These  conditions  are 

independent  of  £. 

Because  of  (15)  we  can  consider  the  change  of  probability  given 
by  the  “Girsanov  transformation 

(«)  #|  ■ 

F 

Under  the  probability  P^  the  processes  w^(t)  and  v^(t)  are  independent 
standard  Wiener  processes.  Since  w  (t)  and  v  (t)  are  independent  of 
Plunder  P  Cis  independent  of  w^(t)  and  v*^(t).  Further,  since 
|l^^(t),  }  is  a  martingale,  C  has  the  same  distribution  under  P  as 

under  P. 


7.3.  The  Filtering  Problem 


In  the  space  (  S}  ,P,P^,P^)  we  can  write 


dx 


(24) 


g(— )  X  dt  +  o(— )dw 
x^(0)  -  C 


dz'^  -  h(^)  x^dt  +  dv^(t) 


e  e  •  t 

where  w  and  v  are  standard  F  -  Viener  processes  which  are  mutually 
independent.  Moreover,  C  is  a  F^  -  Gaussian  random  variable  with  mean 
X  Q  and  covariance  matrix  Pq .  The  filtering  problem  associated  with 

(24)  consists  in  computing 

(25)  ir®^(t)(4»)  -  E^[tj»(x®(t))  I  Z*^l 

for  any  Borel  bounded  test  function  on  R  .  It  is  easy  to  check  that 


ir''(t)(./;) 


(26) 


P^(t)M 

P^(t)(l) 


where 

(27)  P^(t)(«|»)  -  E[i|.(x''(t))u^(t)|  Z*"] 

Our  purpose  here  is  to  study  the  behavior  of  this  quantity  as  E-+-  0. 
In  subsequent  arguments  it  is  useful  to  have  the  bound 

(28)  •  E  <  C. 

To  ensure  this,  we  proceed  as  follows:  For  s  >  1  we  write 
y  (T)  ■  exp {2  /J  (hx^  •  dz  +  (o  ^g)x^  •  dw 

-2s  /J  I  (h+a“^g)x^|^dc}  •  exp  {(2s-l)  /T  |  (h+a‘^g)x'^  |  ^dt } 


(29) 


From  this  ve  have 


(30)  E  u^(T)^  1  (E  exp  /q  l(h+o'-"g)x^l^dt})  s 


(2s-l)  .T  I I  2  , 


Note  that  3(28-l)/(8-l}  has  a  mlnlaun  on  [l,<»)  at  some  ^  >  1.  Thus, 
it  suffices  to  check  that 


s  (2S--1)  ^  _ 

(31)  E  exp{-^^- — :j -  Tl(h+a  ^g)x^(t)|  }<«,  te[0,T] 


This  is  similar  to  (14)  except  that  the  parcuneter  S  is  fixed.  Tcdclng 


s  (2s^-l)  ,  2 

(32)  6  t  -  T  II  h+o  ^g  11^ 


we  require  (22)  which  reads 

1  >  4  11  a  11  T^  11  h+o"^g  11  ^  -  «o 


(1  -  6q)1  >  26Po. 

These  conditions  restrict  the  size  of  T,  and  the  extent  to  which  they 
are  necessary  is  not  clear. 


7.4.  A  Duality  Form  and  an  Expression  for  the  Conditional  Density 


By  introducing  a  certain  duality  formula  it  is  possible  to  obtain 
an  expression  for  the  conditional  density  which  is  convenient  for  the 
homogenization  and  convergence  analysis. 


Let  B  be  a  determiniatic  function  in  L  (0,T;B  )  and 
(34)  P(t)  •  exp  {  /q  6  •  dz  -  %  /J  jel^ds} 

It  is  known  that  VT,  the  set  of  random  variables,  {^(t)},  obtained  by 
varying  B  in  l"(0,T;R*^)  is  dense  in  L  (  n,^,P;S^). 


Let  '1'  be  a  smooth,  bounded  function  on  R  and  let  ^  (t)  be  a 
smooth,  bounded  deterministic  function  on  [0,t]  with  values  in  R  .  Ve 
Introduce  the  deterministic  function  v£(z,t)  which  is  the  solution  of 

Ir  * 

v®(x,T)  -  (x)  ,  T  ^  t  ^  0 


Because  the  Coefficients  are  smooth,  (35)  has  a  solution  in 
X  [0,t]).  Moreover,  it  satisfies  the  growth  conditions 


(36)  jv*^  {x,t)l<  Cg  exp  6|x|^ 

|Dv‘^(x,t)l<  exp  25 1x1^ 

where  5>  Ocan  be  chosen  arbitrarily  small.  Note  that  the  first 

constant  C  ,  in  (36)  can  be  chosen  independent  of  e,  but  not  5. 

0 

One  way  to  verify  (36)  is  to  use  a  probabilistic  formula  for 
v^(x,t).  Consider  the  equation 


(37) 


■  o o -.■  -.^V" o ■  V w."Tiw - -  — 


e 


n 


■ 


1 1 


y. 


p 


i  :■:• 


i.*' 


I  ^ 


«  '•/ 


c 


e  c 

dx^  -  g(~)  x^dt  +  o(^)db 

c  & 


x^(t) 


on  a  probabilitj  space  (not  necessarily  the  original  one)  where  b(s) 
is  a  standard  Wiener  process*  Then 


(58)  v*^(x.t)  -  E  {.|)(x^(T))  exp  /J  h(^)x'^6ds} 

t  ^ 


Therefore , 


Iv'^Cx.t)]  ^  K  E  exp  Clx^(s)|  ds 


(39) 


<  K,  fl  E  exp  6lx^(s)|^  ds 
—  0  t 


where  6>  0  may  be  chosen  arbitrarily  small.  A  calculation  similar  to 
(18)  shows  that 


(40) 


E  exp  6  |x®(t)|^^  kj(0)exp  pj(0)|x|' 


where 


pj(t)  -  6,  t  ^  s 

lUII  +  2||sll  p!  -  0 


ll|_+  2p'  n  lUII  -  0  .  k‘<t)  .  1 


Now 


pj(s) 


exp[-2  ||g||  (t-s)  -  4  ||a||5(t-s)]  (l.-exp[-,2  iL£iJ_Lt-s 

2  l|g||  (t-s) 
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L*Jl3 


(42)  2ll8ll(t-s)  exp[-2  UgH  (t-s)l 
1-exp [-2  ||g||  (t-s)] 


2l|g||  T  exp  t-2lj: 
1-exp  [2  llgll  T] 


li 


Since  the  function  xexp(-x)/[l-exp(-x)]  is  decreasing  on 
has 


'] ,  one 


.  _ 26  II  g||  (t-s) _ 

a-.*p[-2  llsll  (t-,)!)  IIbI  (t-Ol  .4 

l-exp[-2  ||g||  (t-s)] 


If  we  choose  6  >  0  so  that 


(44)  4  Hall  6T  <  All-S-ll  e-Xp[-2  HgH  T]  , 

1-exp [-2  llgll  T) 


then 


(45)  125(8)!  1  - 11  S  II  T _ _ 

2  llgll  T  exp(-2  llgll  T]-4  ||a||  6T(l-exp[-2  ||g||  T]) 


And  from  this  the  first  estimate  in  (56)  follows. 

To  prove  the  second  estimate  in  (36),  one  may  proceed  by 
differentiating  the  expression  (38).  Namely, 
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(46) 


is-  -  E 

3x^ 


3x^  • 
e, 


-  exp  [  hfex^'edsl 

dx.  t  e 


+  \(»(x*^(T))  (exp(  h(^)x^6ds]) 


,I  ,1  *'*lk  c  .  ,  ,  "t  ,  ,  .  , 

•  <T  ‘jt*  35^  <“>  ■*•  ’ 

and  from  (37), 


3x: 


K 


(470 


K 

3x. 


4 

3x^ 

3x^ 

3x . 
1 

1 

3o 

ki 

db 

e 

3Xj 

**1 

il 

)  ds 


(t) 


"ki 


s  >  t  >  0 


It  follows  from  (47)  that 


3x,  , 

(48)  E .  ( I  3;j^(s)  I  )  1  C(1  +  E  /®  |x^(r)|^dr] 

Hence,  —  C  (1  +  (x(  ). 

^  E  (l-g^ - )  1  C  (1  +  lx|) 

i 

and  from  this  one  can  readily  deduce  the  second  estimate  in  (36). 

Using  the  function  v^(x,t),  it  is  possible  to  obtain  a  convenient 
expression  for  p^(t)(4'). 


Proposition  1.  Under  assumptions  (9)  we  have  for  any  t: ,  the  equality 


(50)  E[p'^(T)(i|.)p(T)] 


E[v*=(C.O)] 


where 

(51) 

Proof 

(52) 

But 

(53) 

and 

(54) 

Using 


-  v^(x,0)TTQ(x)dx 

R 

1  T  -1 

“  - ;; - Ir  ®*P  t"**  (x-xJl- 

^  [(2ir)"detPg]'*  u  0  u 

From  (27)  we  have 

Elp‘^(T)(t|.)6(T)]  -  E(4>(x*=(T))u‘'(T)e(T)] 


dv*^(x®(t)  ,t) 


*2  e 

3  V 


‘ij  3x^3Xj 


)dt 


e  1  ®  'v. 

d(u^p)  =  pu  th(  — )  x^  •  dz  +  a  g(  ^)  •  dw] 

C  £ 

c 

+  puB'dz  +  ppB'h(  ^)x^dt 

this  and  (35),  we  have 


d(v^(x'^(t),t)u'^(t)p(t)l  -  u^(t)p(t)[(a*(“)Dv'^  (x'^  (t),t) 


+  v'^Cx^CO.t)  o'^  g(  ^^^)x^t))  •  dw 
+  v^(x^t).t)  (h(^^^)x^(t)  +  6(t))-dz]. 


(55) 


Because  of  the  estimates  (36)  one  can  take  the  expectation  of  the 
stochastic  integrals  obtained  by  integrating  (55  )•  Integrating  and 
taking  the  expectation  gives 

E  y®(e,  0)  -  E(v®(x®(T),T)  u^(X)p(T)1 

which  is  the  desired  result. 

QED. 

Remark.  Note  that  (50)  is  well  defined  if  ¥  is  Borel  bounded  and 
S  e  iT  (0,T;H^).  In  this  case  the  function  is  not  x  [0,t]); 

but  this  is  not  essential  for  the  right  hand  side  of  (50)  to  be  well 
defined.  Thus,  by  regularization,  it  follows  that  (50)  also  holds 
whenT  is  Borel  bounded  and  S  t:  L  (0,T;R^). 

« 

7.5.  Homogenization 

Our  objective  is  to  derive  a  homogenization  representation  of  the 
conditional  distribution  p*tt)(Y)  as  0.  ¥e  shall  begin  by 
considering  the  homogenization  of  (35),  which  is  a  relatively 
classical  problem.  Formally,  the  method  is  as  follows:  Ve  consider 
an  expansion  of  the  form 

(56)  v^(x,t)  -  Vg(x,t)  +  evj^(x,  t)  +  e^v^Cx,  t)  +  v^(x,t) 

Introducing  y  ■  x/ c  and  using  the  expression 


we  obtain 


3v, 


av 


3v^ 


1  •  -2  Hi  +  i£:  +  .„(y) 


— ^  +  — —  +  c  —  + 

at  at  at  at 


Ij'"  3 


+  +  i. 


a^v. 


3^v, 


<y) 


e  “ij'-'''  ay^ayj 


+  2  a.,(y) 


Ij'"  3yj3»j  * 


(58)  +  (y)xj  (e  +  ^)  +  e  ,^hy  (y)«^  6j 


3^, 


+  a.4(y) 


+  2e  a..(y) 


2  ^^^2 

+  E  (y)  jXj^axj 


2  ^^2 

+  8y<y)*i  <'  37-  + 


ay^syj 

^)+  hj^(y)Xj8j 


-A  V 


where  we  have  set 


^2 

.e  ,  »  a  V 

(59)  A  ^ 


3v 


Ij'"  3x^3»j  -  3I;;-  -  '  '’ij<>'>5/l 


with  y  ■  x/  .  We  choose 


(60) 


v^(x,y,t)  -  Vj^(x,t) 


and 


3Vx 


3V 


3v, 


^  +  a,,(y)3^4r  +  8ij(y)>'j 


aT"  3x^3xj 


(61) 


3^. 


+  Vohi.(y)Xj  6^  a.^(y) 


To  deal  with  the  latter,  we  Introduce  m(y)  the  unique  solution  of 


(62)  3^77  ° 


m  periodic  on  Y,  m  >  0,  me  C^,  /ym(y)dy  ■  1  (cf.  Bensouasan, 
Lions,  and  Papanicolaou  1978,  p.  530).  Then  the  solvability 
condition  (Fredholm  Alternative)  for  (61 )  is 

3vo  _  a^v  _  avQ 

It”  ®ij  ax  ax  ®ij  *j  ax.  ° 

i  j  ^ 


Vq(x,T)  -  4'(x),  T  t  ^  0 


where  we  have  set 


(64)  m(y)dy 


and  similarly  defined  g.  .  and  Ti. 

i]  13 


If  we,  in  fact,  choose 


(65)  Vj^(x,t)  2  0 


then  v^(x,t)  is  the  solution  of 


3v  .  . e 


avz 

ay.Sx.  ^ij'^j  ay 7 


-  avo  a  V,  av, 

+  e(  — -  +  a  -  +  gj  .X.  T —  +  V.,  h.  .X.3.) 

'■at  ij  ax^ax^  ^ij  j  ax^  2  ij  j  i 

'Vg 

V  (x,T)  =  0  ,  T  >  t  >  0 


To  estimate  v  ,  we  proceed  as  follows:  First,  we  derive  an 
explicit  formula  for  'i^Cxft)  which  is  similar  to  (38).  Consider  the 
Gaussian  process 

(67)  d?  -  ic  dt  +  odb  5(t)  -  x 
_  _  1 

where  a  ~  (2a) ^  .  Using  this 

T  _ 

(68)  VQ(x,t)  -  E{i|»  (C^^t(T))  exp(  h 

and  we  can  easily  check  that 

lvQ<x,t)l  <_  exp(6lx|S 

(69) 

l°VQ(x,t)|  <  exp(6lx)^ 

for  some  K  ^  and  any 6  >  0.  An  additional  calculation  shows  that 
3^v 

I  ax  '3x  I  -  h 

i  j 

These  estimates  mean  that 

(70)  I  I  <_  exp(6|x|^). 

From  (6l)-(63)  we  can  assert  that 


(l)  9  is  not  the  average  of  9.  This  is  a  slight  abuse  of  notation. 


V  -J»  v»  *>  V  V  ”1^  *.»  ’-*  ">">  *>  "J*  "J*  V»| 


i  ^ 


f*. 


»i 


3v 


(72) 


V2(x.y.t)  -  x^j(y)  ^— +  n,j(y)x.  — 


+  VpC^jCy)  XjB^, 


for  some  smooth,  bounded  functions  X^^j,  ,  and  on  7*  Since  the 
higher  order  derivatives  of  v^  also  satisfy  the  bounds  (69)  -  (71),  we 
can  deduce  from  (66)  and  (72)  that 


(75)  -|^  +  A^v® 


ef*^v®  (x.T)  -  0 


i  I 


P  c 


I  - 


I  C 


r- 

L, 


where 


(74)  (x»^)|  f.  K  exp(6|x|^) 


Again  considering  (37),  we  can  write 


(75)  V®  (x,t) 


eE  {  f®(x®(s),8)(exp[  h  (^)x^*Bdr])ds} 


And,  by  using  arguments  similar  to  those  which  led  to  the  first 
estimate  in  (36),  we  obtain 


(76)  •  (x,t)|  <_  eK^  exp(6|x|^) 

where  5>  0  can  be  chosen  arbitrarily  small.  Combining  this  estimate 
with  the  expression  (72)  for  Vj,  we  have  proved  the  following: 

Proposition  2.  Under  the  assumption  (l)  we  have  the  estimate 

(77)  |v^(x,t)  -  VQ(x,t)|ieK^  exp(6|x|^) 
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where  6^  0  can  be  chosen  arbitrarily  small. 


By  adapting  this  procedure  we  can  provide  a  similar  analysis  of 
the  homogenization  properties  of  the  conditional  distribution  (27)  in 
the  nonlinear  filtering  problem.  This  is  the  main  result  of  this 
section. 


First,  consider  the  "limiting  filtering  problem"  defined  as 
follows;  Let 


dx  “  g  xdt  +  adw 


(78)  dz  -  h  xdt  +  dv 


x(0)  =  z(0)  =  0  ,  0  <  t  <  T 


and  let 


(79)  p''(T)(.|»)  -  E[4.(x(T))v°(T)|  Z^] 


where 


(80)  v°(t) 


exp  {  /q  hx  •  dz  -  *5  /g  I  ^*1  ‘is) 


in  which  z  is  a  standard  Wiener  process.  ((78)  follows  from  a 
Girsanov  transformation  as  used  in  (24)).  In  fact,  we  have  the 
well-known  formula 


(81)  p®(T)(i<()  -  exp(-p(T)) 


>  (y )  exp  [ ->i(y-x  (T) )  (T)  (y-x  (T)  )  1 

(2iT)"/^(det  P(T))^ 


in  which 


e-_ 


.N 


(82)  P(t)  “  *5  -^0  ~  -^0 


and  z(t)  is  the  state  of  the  Kalman  filter 


dx  ■  gx  dt  +  Pii^  (dz  -  hx  dt) 
x(0)  =  x_ 


P  +  Ph^hP  -  0  0^  -  (gP  +  Pg^'i  =  0 

P(0)  =  P. 


As  In  Proposition  1,  ve  can  show  that 
E[p°(T)(t|»)p(T)]  »  E(Vq(C.  0)1 

(84) 

“  VQ(x,0)xQ(x)dx. 

R 


Using  this,  we  can  state  the  following; 


Theorem.  Under  the  assumptions  (9)  and  (3?)  we  have 


(85)  p^T)0;)  - -  P°(T)(.|») 

G  -*•  0 


weakly  in  (  n,Z^,P)  for  every  hounded,  uniformly  continuous 


Proof.  First  note  that  we  can  assume,  without  loss  of  generality, 
that  ip  is  smooth  and  bounded.  Indeed,  let?^  L^(  n,z'^,P),  then 
using  (28) 


_  •  *!•  * 


where  ^  ('J') 


.  P 

P°(T)(1) 

0 

denotes  the  limit  conditional  probability,  and  E  refers  to  the 
probability  on  Q  for  which  z  satisfies  (78).  Therefore,  we  can  assert 
that 

(t)  ->■  E°Tr°  (T)  {f)C^  vi|;,v?^ 

It  would  be  nice  to  prove  stronger  convergence  results,  but  it  must  be 
kept  in  mind  that  the  processes  (l)  themselves  converge  Just  in  law 
and  not  in  a  stronger  sense  (cf.  Bensoussan,  Lions,  and  Papanicolaou 


8.  Open  Problems  and  Further  Work 


The  two  main  issues  which  we  plan  to  develop  in  subsequent  phases 
of  this  research  program  are  the  integration  of  the  modeling  and 
control  methodologies  along  the  lines  initiated  in  sections  5r  ^  and 
7,  and  the  development  of  numerical  software  to  facilitate  the 
application  of  the  integrated  methodology  to  complex  structures.  The 
specific  issues  which  we  intend  to  address  in  Phase  II  of  this 
research  program  are: 

t.  Homogenization  and  asymptotic  analysis  of  control  (including  state 
estimation)  for  large  lattice  type  structures. 

2.  Viener-Hopf  -  spectral  factorization  methods  for  the  control  of 
complex  space  systems  with  hybrid  (lumped  and  distributed) 
structure. 

3.  Development  of  stabilizing  control  strategies  for  nonlinear 
distributed  models,  including  nonlinear  beam  and  plate  models. 

4.  Synthesis  of  a  design  methodology  for  hybrid  nonlinear  structures, 
including  the  nonlinear  differential  geometric  methods  which  have 
been  used  for  finite  dimensional  control  problems  and  the  (linear) 
methods  which  we  have  developed  for  the  treatment  of  distributed 


linear  models. 


1 .  Asymptotic  Analysis  and  Homogenization  of  Control  Problems 
for  Lattice  Type  Structures. 

In  Phase  I  of  this  program  the  use  of  homogenization  in 
connection  vith  control  system  designs  was  demonstrated  in  the 
analysis  of  two  abstract  control  and  filtering  problems.  This 
analysis  established  the  mathematical  feasibility  of  the  technique. 
Previously,  homogenization  had  been  used  only  for  model  reduction,  and 
it  had  not  been  applied  in  a  control  theoretic  setting.  Since  the 
basic  optimization  techniques,  like  the  Bellman  equation,  are 
inherently  nonlinear,  it  was  not  clear  how  the  methodology  could  be 
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used.  The  feasibility  of  the  homogenization  methodology  has  now  been 
demonstrated  in  the  context  of  abstract  control  problems. 

We  have  also  demonstrated  use  of  the  methodology  in  the 
construction  of  simplified  models  for  certain  kinds  dynamical 
phenomena  propagating  on  lattice  structures.  We  have  treated  heat 
conduction  type  problems,  by  exploiting  the  connection  between  such 
problems  and  an  associated  probabilistic  structure,  and  simple  one 
dimensional  lattice  structures  using  purely  analytical  (PDE)  methods 
for  model  simplification.  The  asymptotic  analysis  method  involves  a 
study  of  the  convergence  of  the  resolvents  of  certain  operators  using 
the  theory  of  "correctors"  introduced  for  this  class  of  problems  by 
Bensoussan,  Lions,  and  Papanicolaou  (1978).  The  method  handles  the 
transition  from  the  "discrete"  operators  characterizing  lattice  type 
structures  to  the  PDE  operators  characterizing  the  continuum 
approximations  of  the  structures.  It  also  ^  produces  the  natural 
continuum  model  in  the  course  of  the  asymptotic  analysis;  that  is,  it 


is  not  necessary  to  postulate  the  model  a  priori.  For  ezcunple, 
homogenization  of  the  one  dimensional  lattice  structure  in  section  5 
produced  the  Timoshenko  beam  model  rather  than  the  Euler  beam  model 
which  one  might  have  expected  from  the  symmetry  of  the  orginal 
formulation. 

In  the  second  phase  of  this  program  we  shall  develop  the 
methodology  to  treat  the  combined  problem  of  modeling  and  control  of 
lattice  type  structures;  that  is,  we  shall  develop  the  homogenization 
-  optimization  methodology  described  in  section  6  to  treat  realistic 
models  of  the  dynamical  control  of  large  lattice  and  plate  structures. 

I 

2.  Wiener«Hopf  -  spectral  factorization  methods  for  the  control 
of  complex  space  systems  with  hybrid  (lumped  and  distributed) 
structure. 

Thus  far  we  have  applied  our  control  theory  only  to  purely 
distributed  models  of  a  very  simple  character.  While  it  is  clear  that 
the  methods  can  be  used  for  the  design  and  analyis  of  control  systems 
for  structures  with  both  distributed  and  lumped  components,  it  would 
be  useful  to  treat  a  problem  Including  both  kinds  of  elements  with  an 
overall  linear  model.  The  NASA  challenge  problem  is  of  this  type  and 
we  shall  consider  it,  adapting  our  frequency  domain  methods  as 
required  to  carry  out  the  design.  Since  this  problem  is  being 
considered  by  several  researchers  using  a  variety  of  methods, 
comparison  of  results  and  capabilities  will  be  possible.  While  there 
are  no  conceptual  problems  in  this  extension  of  our  methods,  we  do 
expect  to  encounter  challenging  numerical  problems. 


Development  of  etabilieing  control  strategies  for  nonlinear 
distributed  models,  including  nonlinear  beam  and  plate  models. 

Many  of  the  applications  of  large  space  structures  require 
maneuvers  which  cannot  be  faithfully  described  by  linear  models.  For 
example,  large  attitude  excursions  of  telescope  and  antenna  structures 
require  equations  for  the  evolution  of  the  Euler  angles  through  the 
course  of  the  maneuver.  We  plan  to  extend  our  methods  to  treat 
certain  aspects  of  this  class  of  problems.  It  will  be  necessary  to 
use  differential  geometric  methods  to  describe  the  global  dynamics  of 
the  system  undergoing  large  angle  maneuvers.  Recent  work  by  Baillieul 
(I983)f  B1  Baraka  and  Krishnaprasad  (1984),  among  others  has  led  to  a 
theory  for  the  attitude  dynamics  for  articulated  structures.  Elements 
of  this  theory  in  combination  with  the  methods  for  the  control  of 
distributed  systems  which  have  been  used  in  Phase  I  of  this  project 
should  be  a  useful  starting  point  for  the  development  of  a 
comprehensive  theory  for  large  scale  motions  of  complex,  distributed 
structures. 

Specific  issues  to  be  addressed  include  the  use  of  stabilization 
techniques  for  semilinear  distributed  systems.  These  are  systems  in 
which  the  controls  enter  the  dynamics  by  multipling  the  state.  Common 
examples  include  the  dynamics  of  a  beam  in  which  the  applied  load  can 
be  manipulated.  When  the  load  is  used  as  a  feedback  control,  the 
system  is  nonlinear,  and  the  theory  of  nonlinear  semigroups  is  the 
most  convenient  setting  for  the  analysis.  In  a  series  of  papers  Ball 
and  Slemrod  have  derived  conditions  for  the  stabilization  of  such 
systems.  In  particular,  they  show  that  stabilization  of  the  Euler 


beam  by  semilinear  feedback  is  a  delicate  problem,  and  that  the  most 
natural  conditions  tend  to  lead  to  a  weak  form  of  stability. 


In  the  second  phase  of  this  program  we  plan  to  combine  this 
theory  with  the  corresponding  theory  for  the  stabilization  of  finite 
dimensional  nonlinear  systems  undergoing  large  attitude  motions. 


4.  Synthesis  of  the  Nonlinear  and  Distributed  Design  Methods 


Resolution  of  the  problems  in  1.-3*,  will  lead  to  a  design 
procedure  for  (a  class  of)  control  systems  for  large  space  structures. 
Work  will  be  necessary  to  unify  the  various  methods  into  a  software 
system  for  computer-aided-design.  We  shall  use  software  systems  for 
symbolic  manipulation  (either  HACSYMA  or  SHP)  to  implement  the  complex 
analysis  involved  in  the  initial  reduction  of  the  modeling  and  control 
problem  (for  example,  by  carrying  out  an  asymptotic  analysis  in  the 
context  of  the  control  system  design).  This  will  pennit  us  to  base 
the  selection  of  numerical  routines  for  the  implementation  of  the 
models  and  control  laws  on  simplified  structural  models.  This  will  in 
turn  reduce  the  number  and  diversity  of  costly  computer  runs  which 
must  be  made  with  conventional  design  tools. 
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APPEinSXX  A 


Experinental  Results  For  The  Culcr  Beam 

The  plots  appear  in  groups.  A  first  page  of  each  group  contains  parameters 
of  the  run.  For  example  the  first  group  is  preceeded  by  a  page  vith  the 
following  parameters.  CONTROL:  .2,. 5. >.8.  The  number  designate  locations 
of  point  actuators.  Similar  remarks  apply  for  the  OBSERV,  where  this 
OBSERV  designates  the  points  on  the  beam  whose  displacement  is  penalized  in 
the  cost  criterion.  FREQUBNCY  RANGE  is  the  range  of  frequency  over  which 
we  evaluate  the  transfer  function  matrix,  the  spectral  factors,  and  the 
resolvent.  RELATIVE  WEIGHT  OP  X  VS  U  is  the  weight  of  the  norm  of 
observation  vector,  assuming  the  weight  of  the  control  vector  to  be  one. 
The  name  of  the  file  is  a  working  variable.  WHAT  NODE:  here  we  list  all 
the  special  modes  that  are  used  to  displace  the  beam  (one  at  a  time),  and 
the  next  line  provides  the  amplitude  of  this  disturbance.  FEEDBACK  GAIN: 
is  -1  for  the  optimal  gain  and  it  is  0  for  the  open  loop.  The  meaning  of 
the  plot  titles  is  as  follows:  "Beam  at  x,  yth  mode,"  means  that  the  beam 
is  initially  displaced  by  the  y-th  special  mode,  and  the  deflection  of  the 
beam  is  observed  as  a  function  of  time  at  point  x.  On  the  control  plot  we 
indicate  the  position  of  the  point  control  and  the  time  evolution  of  the 
control  at  that  point. 
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